added uniformly distributed points calculation for PCS (also usefull for BTA & Z-1000)

This commit is contained in:
Edward Emelianov 2022-10-05 17:32:35 +03:00
parent b26a4fa173
commit c520f1146c
13 changed files with 592 additions and 0 deletions

57
calc_points/Makefile Normal file
View File

@ -0,0 +1,57 @@
# run `make DEF=...` to add extra defines
PROGRAM := calcpoints
LDFLAGS := -fdata-sections -ffunction-sections -Wl,--gc-sections -Wl,--discard-all
LDFLAGS += -lusefull_macros -lm
SRCS := $(wildcard *.c)
DEFINES := $(DEF) -D_GNU_SOURCE -D_XOPEN_SOURCE=1111
OBJDIR := mk
CFLAGS += -O2 -Wall -Wextra -Wno-trampolines -std=gnu99
OBJS := $(addprefix $(OBJDIR)/, $(SRCS:%.c=%.o))
DEPS := $(OBJS:.o=.d)
TARGFILE := $(OBJDIR)/TARGET
CC = gcc
#TARGET := RELEASE
ifeq ($(shell test -e $(TARGFILE) && echo -n yes),yes)
TARGET := $(file < $(TARGFILE))
else
TARGET := RELEASE
endif
ifeq ($(TARGET), DEBUG)
.DEFAULT_GOAL := debug
endif
release: $(PROGRAM)
debug: CFLAGS += -DEBUG -Werror
debug: TARGET := DEBUG
debug: $(PROGRAM)
$(TARGFILE): $(OBJDIR)
@echo -e "\t\tTARGET: $(TARGET)"
@echo "$(TARGET)" > $(TARGFILE)
$(PROGRAM) : $(TARGFILE) $(OBJS)
@echo -e "\t\tLD $(PROGRAM)"
$(CC) $(LDFLAGS) $(OBJS) -o $(PROGRAM)
$(OBJDIR):
@mkdir $(OBJDIR)
ifneq ($(MAKECMDGOALS),clean)
-include $(DEPS)
endif
$(OBJDIR)/%.o: %.c
@echo -e "\t\tCC $<"
$(CC) -MD -c $(LDFLAGS) $(CFLAGS) $(DEFINES) -o $@ $<
clean:
@echo -e "\t\tCLEAN"
@rm -rf $(OBJDIR) 2>/dev/null || true
xclean: clean
@rm -f $(PROGRAM)
.PHONY: clean xclean

19
calc_points/Readme Normal file
View File

@ -0,0 +1,19 @@
Calculate coordinates of points equally distributed by hemisphere.
Usage: calcpoints [args]
Where args are:
-1, --hideA hide first column (A)
-2, --hideZ hide second column (Z)
-3, --hideHA hide third column (HA)
-4, --hideDec hide fourth column (Dec)
-Z, --maxz=arg maximal Z (degrees) (default: 75.)
-c, --scoord=arg sort by this coordinate (A, Z, HA, Dec) (default: Z)
-d, --delimeter=arg coordinates delimeter string (default: ':')
-h, --help show this help
-n, --npts=arg max amount of points (default: 100)
-o, --output=arg output file name
-s, --sorting=arg sorting order (none, positive, negative)
-z, --minz=arg minimal Z (degrees) (default: 0.)

86
calc_points/cmdlnopts.c Normal file
View File

@ -0,0 +1,86 @@
/*
* This file is part of the uniformdistr project.
* Copyright 2022 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>
#include <stdio.h>
#include <string.h>
#include <strings.h>
#include <math.h>
#include "cmdlnopts.h"
#include "usefull_macros.h"
/*
* here are global parameters initialisation
*/
static int help;
// default global parameters
static glob_pars G = {
.delimeter = ":",
.Npts = 100,
.Zmax = 75.,
.sortcoord = "Z",
};
/*
* Define command line options by filling structure:
* name has_arg flag val type argptr help
*/
static myoption cmdlnopts[] = {
// common options
{"help", NO_ARGS, NULL, 'h', arg_int, APTR(&help), "show this help"},
{"delimeter",NEED_ARG, NULL, 'd', arg_string, APTR(&G.delimeter), "coordinates delimeter string (default: ':')"},
{"output", NEED_ARG, NULL, 'o', arg_string, APTR(&G.outfile), "output file name"},
{"npts", NEED_ARG, NULL, 'n', arg_int, APTR(&G.Npts), "max amount of points (default: 100)"},
{"minz", NEED_ARG, NULL, 'z', arg_double, APTR(&G.Zmin), "minimal Z (degrees) (default: 0.)"},
{"maxz", NEED_ARG, NULL, 'Z', arg_double, APTR(&G.Zmax), "maximal Z (degrees) (default: 75.)"},
{"sorting", NEED_ARG, NULL, 's', arg_string, APTR(&G.sorting), "sorting order (none, positive, negative)"},
{"scoord", NEED_ARG, NULL, 'c', arg_string, APTR(&G.sortcoord), "sort by this coordinate (A, Z, HA, Dec) (default: Z)"},
{"hideA", NO_ARGS, NULL, '1', arg_none, APTR(&G.hide[0]), "hide first column (A)"},
{"hideZ", NO_ARGS, NULL, '2', arg_none, APTR(&G.hide[1]), "hide second column (Z)"},
{"hideHA", NO_ARGS, NULL, '3', arg_none, APTR(&G.hide[2]), "hide third column (HA)"},
{"hideDec", NO_ARGS, NULL, '4', arg_none, APTR(&G.hide[3]), "hide fourth column (Dec)"},
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){
int i;
size_t hlen = 1024;
char helpstring[1024], *hptr = helpstring;
snprintf(hptr, hlen, "Usage: %%s [args]\n\n\tWhere args are:\n");
// format of help: "Usage: progname [args]\n"
change_helpstring(helpstring);
// parse arguments
parseargs(&argc, &argv, cmdlnopts);
if(help) showhelp(-1, cmdlnopts);
if(argc > 0){
G.rest_pars_num = argc;
G.rest_pars = MALLOC(char *, argc);
for(i = 0; i < argc; i++)
G.rest_pars[i] = strdup(argv[i]);
}
return &G;
}

39
calc_points/cmdlnopts.h Normal file
View File

@ -0,0 +1,39 @@
/*
* This file is part of the uniformdistr project.
* Copyright 2022 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
/*
* here are some typedef's for global data
*/
typedef struct{
char *outfile; // output file name
char *delimeter; // coordinates delimeter: HH:MM:SS.S, HH/MM/SS.S etc
char *sorting; // soring function (none, increasing, decreasing)
char *sortcoord; // coordinate to sort by
int rest_pars_num; // number of rest parameters
int Npts; // max amount of points
int hide[4]; // hide n-th column
double Zmin; // minZ (degr)
double Zmax; // maxZ (degr)
char** rest_pars; // the rest parameters: array of char*
} glob_pars;
glob_pars *parse_args(int argc, char **argv);

128
calc_points/main.c Normal file
View File

@ -0,0 +1,128 @@
/*
* This file is part of the uniformdistr project.
* Copyright 2022 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 <stdarg.h>
#include <stdio.h>
#include <usefull_macros.h>
#include "cmdlnopts.h"
#include "uniform.h"
// return sign of angle (+1 - positive, -1 - negative, 0 - error)
static int deg2dms(double deg, int *d, int *m, int *s){
if(!d || !m || !s) return 0;
int sign = 1;
if(deg < 0.){
sign = -1;
deg = -deg;
}
// add 0.5'' for right roundness
deg += 0.5/3600.;
*d = (int) deg;
deg = (deg - *d) * 60.; // minutes
*m = (int) deg;
deg = (deg - *m) * 60.;
*s = (int)deg;
return sign;
}
/*
// for hour angle in hour format
static int deg2hms(double deg, int *d, int *m, int *s){
// if(deg < 0.) deg += 360.;
deg /= 15.;
return deg2dms(deg, d, m, s);
}*/
static void tee(FILE *f, const char *fmt, ...){
va_list ar;
va_start(ar, fmt);
vprintf(fmt, ar);
va_end(ar);
if(f){
va_start(ar, fmt);
vfprintf(f, fmt, ar);
va_end(ar);
}
}
/**
* @brief savepoints - print points and save to file (if f!=NULL)
* @param f - file
* @param pts - points
* @param N - amount
* @param delim - angle delimeter string
* @param mask - bit mask of columns (e.g. 1100 - show A/Z and hile HA/Dec)
*/
static void savepoints(FILE *f, point *pts, int N, char *delim, int mask){
char str[4][32];
if(!pts || !delim || N < 1) return;
const char *headers[4] = {"A", "Z", "HA", "Dec"};
int show[4] = {0};
for(int i = 0; i < 4; ++i) if(mask & (1<<i)) show[i] = 1;
tee(f, "%-15s", "#");
for(int i = 0; i < 4; ++i){
// check which columns we need
if(show[i]) tee(f, "%-14s", headers[i]);
}
tee(f, "\n");
for(int i = 0; i < N; ++i, ++pts){
char *sgn[4] = {"", "", "", ""};
int d[4], m[4], s[4];
if(deg2dms(pts->A, &d[0], &m[0], &s[0]) < 0) sgn[0] = "-";
if(deg2dms(pts->Z, &d[1], &m[1], &s[1]) < 0) sgn[1] = "-";
if(deg2dms(pts->HA/15., &d[2], &m[2], &s[2]) < 0) sgn[2] = "-";
if(deg2dms(pts->Dec, &d[3], &m[3], &s[3]) < 0) sgn[3] = "-";
for(int i = 0; i < 4; ++i)
if(show[i]) snprintf(str[i], 31, "%s%02d%s%02d%s%02d",
sgn[i], d[i], delim, m[i], delim, s[i]);
tee(f, "%-6d", i);
for(int i = 0; i < 4; ++i){
if(show[i]) tee(f, "%14s", str[i]);
}
tee(f, "\n");
}
}
int main(int argc, char **argv){
initial_setup();
glob_pars *G = parse_args(argc, argv);
FILE *f = NULL;
if(G->outfile){
f = fopen(G->outfile, "w");
if(!f) ERR("Can't open %s", G->outfile);
}
if(G->Npts < 10) ERRX("Need at least 10 points");
if(!set_Zlimits(G->Zmin, G->Zmax)) ERRX("Wrong Z limits");
if(G->sorting && !set_sorting(G->sorting, G->sortcoord)){
show_sorting_help();
return 1;
}
int mask = 0x0f;
for(int i = 0; i < 4; ++i){
if(G->hide[i]) mask &= ~(1<<i);
}
if(mask == 0) ERRX("You can't hide ALL columns");
int N = G->Npts;
point *points = getPoints(&N);
if(!points) ERRX("Can't calculate");
green("%d -> %d\n", G->Npts, N);
savepoints(f, points, N, G->delimeter, mask);
if(f) fclose(f);
return 0;
}

212
calc_points/uniform.c Normal file
View File

@ -0,0 +1,212 @@
/*
* This file is part of the uniformdistr project.
* Copyright 2022 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 <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <usefull_macros.h>
#include "uniform.h"
static double Zmin, Zmax;
typedef enum{
SORT_NONE,
SORT_POS,
SORT_NEG,
SORT_AMOUNT
} t_sorting;
typedef struct{
const char *alg;
const char *help;
} strpair;
static const strpair sorthelp[SORT_AMOUNT] = {
[SORT_NONE] = {"none", "don't sort"},
[SORT_POS] = {"positive", "sort in increasing order"},
[SORT_NEG] = {"negative", "sort in decreasing order"},
};
typedef enum{
COORD_A,
COORD_Z,
COORD_HA,
COORD_DEC,
COORD_AMOUNT
} t_coord;
static const char *coordname[COORD_AMOUNT] = {
[COORD_A] = "A",
[COORD_Z] = "Z",
[COORD_HA] = "HA",
[COORD_DEC] = "DEC"
};
// convert angle in radians into 0..360 degrees
static double deg360(double rad){
rad *= 180. / M_PI;
int n = (int)(rad / 360.);
rad -= 360. * n;
if(rad < 0.) rad += 360.;
return rad;
}
// convert angle in radians into -180..180 degrees
static double deg180(double rad){
rad = deg360(rad);
if(rad > 180.) rad -= 360.;
return rad;
}
// compare coordinates - for sorting
static int spos(double a1, double a2){
if(a1 > a2) return 1;
else return -1;
}
static int sneg(double a1, double a2){
if(a1 < a2) return 1;
else return -1;
}
// default compare functions
static int (*compf)(const void *, const void *) = NULL;
// compare points - for sorting
static int compapos(const void *a1, const void *a2){
return spos(((point*)a1)->A, ((point*)a2)->A);
}
static int companeg(const void *a1, const void *a2){
return sneg(((point*)a1)->A, ((point*)a2)->A);
}
static int compzpos(const void *a1, const void *a2){
return spos(((point*)a1)->Z, ((point*)a2)->Z);
}
static int compzneg(const void *a1, const void *a2){
return sneg(((point*)a1)->Z, ((point*)a2)->Z);
}
static int comphapos(const void *a1, const void *a2){
return spos(((point*)a1)->HA, ((point*)a2)->HA);
}
static int comphaneg(const void *a1, const void *a2){
return sneg(((point*)a1)->HA, ((point*)a2)->HA);
}
static int compdecpos(const void *a1, const void *a2){
return spos(((point*)a1)->Dec, ((point*)a2)->Dec);
}
static int compdecneg(const void *a1, const void *a2){
return sneg(((point*)a1)->Dec, ((point*)a2)->Dec);
}
static int (*comparray[COORD_AMOUNT][SORT_AMOUNT])(const void *, const void *) = {
[COORD_A] = {[SORT_NONE] = NULL, [SORT_POS] = compapos, [SORT_NEG] = companeg},
[COORD_Z] = {[SORT_NONE] = NULL, [SORT_POS] = compzpos, [SORT_NEG] = compzneg},
[COORD_HA] = {[SORT_NONE] = NULL, [SORT_POS] = comphapos, [SORT_NEG] = comphaneg},
[COORD_DEC] = {[SORT_NONE] = NULL, [SORT_POS] = compdecpos, [SORT_NEG] = compdecneg},
};
// convert horizontal to equatorial
static void hor2eq(point *p){
double alt_s, alt_c, A_s, A_c; // sin/cos of alt and Az
sincos((90. - p->Z)*M_PI/180., &alt_s, &alt_c);
sincos(p->A*M_PI/180., &A_s, &A_c);
double sind = SIN_LAT_OBS * alt_s + COS_LAT_OBS * alt_c * A_c;
double d = asin(sind), cosd = cos(d);
p->Dec = d *180./M_PI;
double x = (alt_s - sind * SIN_LAT_OBS) / (cosd * COS_LAT_OBS);
double y = -A_s * alt_c / cosd;
p->HA = atan2(y, x) *180./M_PI;
if(p->HA < 0.) p->HA += 360.;
}
/**
* @brief getHpoints - calculate equally distributed points A/Z (degrees)
* @param N (io) - max number of points (i), actual number of points (o)
* @return array of unsorted coordinates
*/
static point* getHpoints(int *N){
if(!N || *N < 1) return NULL;
int n = *N, total = 0;
point* points = MALLOC(point, n), *pptr = points;
double ang = M_PI * (1. + sqrt(5));
for(int i = 0; i < n; ++i){
// we need only one hemisphere, so instead of acos(1. - 2.*(i+0.5)/n) get
double phi = acos(1. - (i+0.5)/n);
double theta = ang * i;
//DBG("phi = %g (%g), theta = %g (%g)", phi, deg360(phi), theta, deg360(theta));
if(phi > M_PI_2) break;
double Z = deg360(phi);
if(Z < Zmin || Z > Zmax) continue;
++total;
pptr->A = deg180(theta);
pptr->Z = Z;
//BG("%d:\t%g / %g", total, pptr->A, pptr->Z);
++pptr;
}
*N = total;
return points;
}
int set_Zlimits(double minz, double maxz){
if(minz > 80. || minz < 0.) return FALSE;
if(maxz > 90. || maxz < 10.) return FALSE;
Zmin = minz;
Zmax = maxz;
return TRUE;
}
point *getPoints(int *N){
point *points = getHpoints(N);
if(!points) return NULL;
// convert A,Z -> Ha,Dec (if sorting by HA/Dec?)
for(int i = 0; i < *N; ++i)
hor2eq(&points[i]);
if(compf){ // sort data
DBG("sort");
qsort(points, *N, sizeof(point), compf);
}
return points;
}
// set sorting parameters
int set_sorting(char *param, char *coord){
DBG("Try to set sorting '%s' by '%s'", param, coord);
int x = 0, c = 0;
for(; c < COORD_AMOUNT; ++c){
if(strcasecmp(coordname[c], coord) == 0) break;
}
if(c == COORD_AMOUNT) return FALSE;
for(; x < SORT_AMOUNT; ++x){
//DBG("x=%d, cmp %s %s", x, sorthelp[x].alg, param);
if(strncmp(param, sorthelp[x].alg, strlen(param)) == 0) break;
}
if(x == SORT_AMOUNT) return FALSE;
DBG("x=%d, c=%d", x,c);
compf = comparray[c][x];
return TRUE;
}
void show_sorting_help(){
fprintf(stderr, "Sorting algorythms:\n");
for(int i = 0; i < SORT_AMOUNT; ++i){
fprintf(stderr, "\t%s - %s\n", sorthelp[i].alg, sorthelp[i].help);
}
fprintf(stderr, "Sorting coordinates:\n");
for(int i = 0; i < COORD_AMOUNT; ++i){
fprintf(stderr, "\t%s\n", coordname[i]);
}
}

39
calc_points/uniform.h Normal file
View File

@ -0,0 +1,39 @@
/*
* This file is part of the uniformdistr project.
* Copyright 2022 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
// SAO coordinates: longitude 41degr 26' 29.175'', latitude 43degr 39' 12.7''
#define LONG_OBS 41.44143375
#define LAT_OBS 43.6535278
#define COS_LAT_OBS 0.723527278
#define SIN_LAT_OBS 0.690295790
typedef struct{
double A;
double Z;
double HA;
double Dec;
} point;
point *getPoints(int *N);
int set_sorting(char *param, char *coord);
int set_Zlimits(double minz, double maxz);
void show_sorting_help();

View File

@ -0,0 +1 @@
-std=c17

View File

@ -0,0 +1,3 @@
#define _GNU_SOURCE
#define _XOPEN_SOURCE 1111

View File

@ -0,0 +1 @@
[General]

View File

@ -0,0 +1 @@
-std=c++17

View File

@ -0,0 +1,5 @@
cmdlnopts.c
cmdlnopts.h
main.c
uniform.c
uniform.h

View File

@ -0,0 +1 @@
.