copy

Makefile

@@ -0,0 +1,24 @@
+# run `make DEF=...` to add extra defines
+OUT = zerndeco
+LDFLAGS = -lm -lgsl -lgslcblas -lpng -lcfitsio -fopenmp
+SRCS = $(wildcard *.c)
+CC = gcc
+DEFINES = $(DEF) -D_GNU_SOURCE -DEBUG
+CFLAGS = -std=gnu99 -fopenmp -Wall -Werror -Wextra $(DEFINES)
+OBJS = $(SRCS:.c=.o)
+
+all : $(OUT)
+
+$(OUT) : $(OBJS)
+	$(CC) $(OBJS) $(LDFLAGS) -o $(OUT)
+
+$(SRCS) : %.c : %.h $(INDEPENDENT_HEADERS)
+	touch $@
+
+%.h: ;
+
+clean:
+	/bin/rm -f *.o *~
+
+depend:
+	$(CC) -MM $(SRCS)

Readme.md

@@ -0,0 +1,9 @@
+Wavefront decomposition by Zernike polynomials
+==============================================
+
+
+Coordinate system - Cartesian (X from left to right, Y from bottom to top)
+On images coordinate system the same, so first byte of image is its left bottom corner!
+
+make:
+DEF=DWF_EPSILON=eps -  set minimal value of vawefront that would mean zero to `eps` (default: 1e-12)

benchmark.c

@@ -0,0 +1,307 @@
+/*
+ * benchmark.c - test of different methods
+ *
+ * Copyright 2016 Edward V. Emelianov <eddy@sao.ru, 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 2 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+#include <stdio.h>
+#include <math.h>
+#include <time.h> // for srand48(time(NULL))
+#include "benchmark.h"
+#include "usefull_macros.h"
+#include "cmdlnopts.h"
+
+void free_wavefront(wavefront **F){
+	FREE((*F)->zdata);
+	FREE((*F)->coordinates);
+	FREE(*F);
+}
+static void free_stat(WF_stat *st){
+	if(!st) return;
+	FREE(st->ZC_sum);
+	FREE(st->ZC_sum2);
+	FREE(st->max_dZC);
+}
+
+void free_bench(WF_benchmark **bench){
+	if(!bench || !*bench) return;
+	WF_benchmark *b = *bench;
+	for(int i = 0; i < WF_bench_size; ++i)
+		free_stat(&(b->stat[i]));
+	FREE(*bench);
+}
+
+/**
+ * calculate Z-coordinate of surface with rectangle nodes by given Zernike coefficients
+ * @param sz     - image size
+ * @param Zidxs  - Zernike coefficients
+ * @param znum   - amount of Zidx
+ */
+double *calc_surface(int sz, double *Zidxs, int znum){
+	if(!Zidxs || znum < 1) return NULL;
+	return Zcompose(znum, Zidxs, sz, sz);
+}
+
+/**
+ * calculate surface with rectangle nodes by given Zernike coefficients using `R`-functions
+ * @param sz     - image size
+ * @param Zidxs  - Zernike coefficients
+ * @param znum   - amount of Zidx
+ */
+wavefront *calc_surfaceR(int sz, double *Zidxs, int znum){
+	int G2 = sz * sz;
+	if(!Zidxs || znum < 1) return NULL;
+	wavefront *F = MALLOC(wavefront, 1);
+	F->size = G2;
+	F->coordinates = MALLOC(polar, G2);
+	polar *rptr = F->coordinates;
+	double Stp = 2./((double)sz - 1.);
+	for(int j = 0; j < sz; ++j){
+		double y = ((double) j) * Stp - 1.;
+		for(int i = 0; i < sz; ++i, ++rptr){
+			double x = ((double) i)* Stp - 1.;
+			double R2 = x*x + y*y;
+			rptr->r = sqrt(R2);
+			rptr->theta = atan2(y, x);
+		}
+	}
+	F->zdata = ZcomposeR(znum, Zidxs, G2, F->coordinates);
+	return F;
+}
+
+/**
+ * calculate surface with scattered nodes by given Zernike coefficients using `R`-functions
+ * @param sz     - amount of points
+ * @param points - given nodes coordinates
+ * @param Zidxs  - Zernike coefficients
+ * @param znum   - amount of Zidx
+ */
+wavefront *calc_surfaceRS(int sz, polar *points, double *Zidxs, int znum){
+	if(!Zidxs || !points || znum < 1 || sz < 1) return NULL;
+	wavefront *F = MALLOC(wavefront, 1);
+	F->size = sz;
+	F->coordinates = MALLOC(polar, sz);
+	memcpy(F->coordinates, points, sizeof(polar)*sz);
+	F->zdata = ZcomposeR(znum, Zidxs, sz, F->coordinates);
+	return F;
+}
+
+/**
+ * calculate coordinates of wavefront for BTA Hartmann mask
+ */
+polar *calc_BTA_Hpoints(int *sz){
+	// radius on Hartmann mask (in mm) - for normalize coordinates @ [-1, 1]
+	double rmax = MIR_R * (1. - HARTMANN_Z/FOCAL_R);
+	double r[] = {175., 247., 295., 340., 379., 414., 448., 478.}; // radii of rings
+	polar *pts = MALLOC(polar, 258), *pt = pts;
+	const double ray_step = M_PI/16., theta0 = M_PI/32; // angle between rays & first ray theta coordinate
+	for(int ring = 0; ring < 8; ++ring){ // main rings
+		for(int ray = 0; ray < 32; ++ray, ++pt){
+			pt->r = r[ring] / rmax;
+			pt->theta = ray * ray_step + theta0;
+		}
+	}
+	// and markers: ring 2 ray 24.5 and ring 7 ray 29.5
+	pt->r = r[2] / rmax; pt->theta = 24.5 * ray_step + theta0;
+	++pt;
+	pt->r = r[7] / rmax; pt->theta = 29.5 * ray_step + theta0;
+	if(sz) *sz = 256;
+	return pts;
+}
+
+static char *bench_names[] = {
+	N_("Direct decomposition (ZdecomposeR)"),
+	N_("Less square decomposition (LS_decompose)"),
+	N_("QR-decomposition (QR_decompose)"),
+	N_("Direct gradient decomposition (directGradZdecomposeR)"),
+	N_("LS gradient decomposition (LS_gradZdecomposeR)")
+};
+
+/**
+ * calculate errors and fill WF_stat fields for single measurement
+ * @param stat    - output stat structure
+ * @param Zidxs   - calculated ZC
+ * @param Zidxs0  - original ZC
+ * @param znum    - size of 'Zidxs'
+ * @param WF0     - original wavefront
+ * @param verbose - ==1 for messages in stdout
+ */
+static void fill_stat(WF_benchmark *bench, WF_bench_type type, double *Zidxs, double *Zidxs0, wavefront *WF0, int verbose){
+	double *zdata0 = WF0->zdata;
+	polar *P = WF0->coordinates;
+	size_t i, sz = WF0->size;
+	WF_stat *stat = &bench->stat[type];
+	int znum = bench->Znum;
+	// test errors by wavefront
+	double *zdata = ZcomposeR(znum, Zidxs, sz, P);
+	double maxd = 0., sum = 0., sum2 = 0.;
+	for(i = 0; i < sz; ++i){
+		double z = zdata0[i];
+		double diff = z - zdata[i];
+		sum += diff;
+		sum2 += diff*diff;
+		if(fabs(z) > WF_EPSILON){
+			diff = fabs(diff / z);
+			if(maxd < diff) maxd = diff;
+		}
+	}
+	stat->max_dWF = maxd;
+	sum /= sz; sum2 /= sz;
+	stat->WF_std = sqrt(sum2 - sum*sum);
+	stat->ZC_sum = MALLOC(double, znum);
+	stat->ZC_sum2 = MALLOC(double, znum);
+	stat->max_dZC = MALLOC(double, znum);
+	if(verbose){
+		printf("%s\n", bench_names[type]);
+		printf("Wavefront error: STD=%g, |max err|=%.3f%%; Zernike:\n#\t   val0   \t    val   \tdiff%%\n", stat->WF_std, maxd);
+	}
+	for(int p = 0; p < znum; ++p){
+		double diff = Zidxs0[p] - Zidxs[p];
+		stat->ZC_sum[p] = diff;
+		stat->ZC_sum2[p] = diff*diff;
+		if(fabs(Zidxs0[p]) > Z_prec) stat->max_dZC[p] = fabs(diff / Zidxs0[p]);
+		if(verbose){
+			printf("%d\t%10f\t%10f\t", p, Zidxs0[p], Zidxs[p]);
+			if(fabs(Zidxs0[p]) > Z_prec) printf("%4.1f\n", stat->max_dZC[p]*100.);
+			else printf("  -  \n");
+		}
+	}
+}
+
+/**
+ * Process a benchmark step with different methods of wavefront restoration
+ * @param Zidxs   - Zernike coefficients
+ * @param znum    - length of Zidxs
+ * @param verbose - ==1 for detailed information on stdout
+ */
+static WF_benchmark *benchmark_step(double *Zidxs0, int znum, int verbose){
+	#define MESG(...) do{if(verbose){printf(__VA_ARGS__); printf("\n");}}while(0)
+	int sz; //G2 = G.imagesize * G.imagesize;
+	wavefront *WF;
+	polar *P;
+	double *Zidxs, *zdata;
+	int N, M, Zsz;
+	convert_Zidx(znum, &N, &M); // calculate max power of Zc
+	WF_benchmark *bench = MALLOC(WF_benchmark, 1);
+	bench->Znum = znum;
+	MESG(_("Direct decomposition of given wavefront"));
+	MESG(_("1. Scattered points for BTA hartmannogram"));
+	P = calc_BTA_Hpoints(&sz);
+	WF = calc_surfaceRS(sz, P, Zidxs0, znum);
+	FREE(P);
+	zdata = WF->zdata;
+	// ZdecomposeR
+	Zidxs = ZdecomposeR(N, sz, WF->coordinates, zdata, &Zsz, NULL);
+	fill_stat(bench, Scatter_direct, Zidxs, Zidxs0, WF, verbose);
+	FREE(Zidxs);
+	// LS_decompose
+	Zidxs = LS_decompose(N, sz, WF->coordinates, zdata, &Zsz, NULL);
+	fill_stat(bench, Scatter_LS, Zidxs, Zidxs0, WF, verbose);
+	FREE(Zidxs);
+	// QR_decompose
+	Zidxs = QR_decompose(N, sz, WF->coordinates, zdata, &Zsz, NULL);
+	fill_stat(bench, Scatter_QR, Zidxs, Zidxs0, WF, verbose);
+	FREE(Zidxs);
+	// directGradZdecomposeR
+	Zidxs0[0] = 0.; // decomposition by gradients cannot know anything about zero coefficient
+	point *grads = directGradZcomposeR(znum, Zidxs0, sz, WF->coordinates);
+	Zidxs = directGradZdecomposeR(N, sz, WF->coordinates, grads, &Zsz, NULL);
+	fill_stat(bench, Scatter_grad, Zidxs, Zidxs0, WF, verbose);
+	FREE(Zidxs);
+	free_wavefront(&WF);
+	/*WF = calc_surface(G.imagesize, Zidxs0, znum);
+	P = WF->coordinates;
+	zdata = WF->zdata;
+	grads = gradZcomposeR(znum, Zidxs0, G2, P);
+	Zidxs = ZdecomposeR(N, G2, WF->coordinates;, zdata, &Zsz, &lastIdx);
+	FREE(grads);
+	*/
+	#undef MESG
+	return bench;
+}
+
+/**
+ * Update values stored in `total` by values of `current`
+ */
+static void update_stat_values(WF_benchmark *total, WF_benchmark *current){
+	++total->measurements;
+	int znum = total->Znum;
+	for(int i = 0; i < WF_bench_size; ++i){
+		WF_stat *out = &total->stat[i], *in = &current->stat[i];
+		if(!out->ZC_sum || !in->ZC_sum) continue; // given field is empty - didn't run benchmark
+		out->WF_std += in->WF_std;
+		if(in->max_dWF > out->max_dWF)
+			out->max_dWF = in->max_dWF;
+		for(int j = 0; j < znum; ++j){
+			out->ZC_sum[j] += in->ZC_sum[j];
+			out->ZC_sum2[j] += in->ZC_sum2[j];
+			if(out->max_dZC[i] < in->max_dZC[i])
+				out->max_dZC[i] = in->max_dZC[i];
+		}
+	}
+}
+
+/**
+ * Process a benchmark with different methods of wavefront restoration
+ * @param Zidxs - user defined Zernike coefficients or NULL
+ * @param znum  - length of Zidxs
+ */
+void do_benchmark(double *Zidxs0, int znum0){
+	int znum = G.bench_maxP;
+	double *Zidxs = MALLOC(double, znum);
+	srand48(time(NULL));
+	inline void gen_zidxs(){ // generate random Z coefficients in interval (-10..10)/p
+		for(int i = 0; i <  znum; ++i)
+			Zidxs[i] = (20.*drand48() - 10.) / ((double)i+1);
+	}
+	WF_benchmark *total = NULL, *current = NULL;
+	if(Zidxs0){ // do only one iteration - for given coefficients
+		current = benchmark_step(Zidxs0, znum0, 1);
+		free_bench(&current);
+		return;
+	}
+	if(G.bench_iter < 2) ERRX(_("Need at least two iterations for benchmark"));
+	gen_zidxs();
+	total = benchmark_step(Zidxs, znum, 0);
+	total->measurements = 1;
+	for(int i = 1; i < G.bench_iter; ++i){
+		gen_zidxs();
+		current = benchmark_step(Zidxs, znum, 0);
+		update_stat_values(total, current);
+		free_bench(&current);
+	}
+	// print stat values
+	double Nm = total->measurements;
+	for(int i = 0; i < WF_bench_size; ++i){
+		WF_stat *stat = &total->stat[i];
+		if(!stat->ZC_sum) continue;
+		printf("%s:\n<STD> = %g\n", bench_names[i], stat->WF_std/Nm);
+		printf("max_dWF=%g%%\n", stat->max_dWF);
+		printf("     p(n,m)\t   Z_std   \tmax_dZC(%%) \tsum2\tsum\n");
+		for(int i = 0; i < znum; ++i){
+			double sum2 = stat->ZC_sum2[i] / Nm;
+			double sum = stat->ZC_sum[i] / Nm;
+			int N,M;
+			convert_Zidx(i, &N, &M);
+			printf("%3d (%2d, %2d)\t%10g\t%10g\t%10g\t%10g\n", i, N, M, sqrt(sum2 - sum*sum),
+				stat->max_dZC[i]*100., sum2, sum);
+		}
+	}
+	free_bench(&total);
+	FREE(Zidxs);
+}
+

benchmark.h

@@ -0,0 +1,71 @@
+/*
+ * benchmark.h
+ *
+ * Copyright 2016 Edward V. Emelianov <eddy@sao.ru, 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 2 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+#pragma once
+#ifndef __BENCHMARK_H__
+#define __BENCHMARK_H__
+
+#include "zernike.h"
+#include "usefull_macros.h"
+
+#ifndef WF_EPSILON
+#define WF_EPSILON (1e-12)
+#endif
+
+typedef struct{
+	size_t size;        // size of arrays
+	double *zdata;      // Z-coordinate of wavefront
+	polar *coordinates; // coordinates on WF (-1..1)
+} wavefront;
+
+// errors of wavefront restoration method
+typedef struct{
+	double WF_std;   // standard deviation for restored wavefront weighted by wavefront values
+	double max_dWF;  // max{|WF - WF_0|/WF_0} if WF_0 > WF_EPSILON
+	double *ZC_sum;  // data for std calculation: sum of differences
+	double *ZC_sum2; //   and sum of their squares
+	double *max_dZC; // max|ZC/ZC_0 - 1|
+} WF_stat;
+
+typedef enum{
+	 Scatter_direct = 0      // ZdecomposeR
+	,Scatter_LS              // LS_decompose
+	,Scatter_QR              // QR_decompose
+	,Scatter_grad            // directGradZdecomposeR
+	,Scatter_LSgrad          // LS_gradZdecomposeR
+	,WF_bench_size
+} WF_bench_type;
+
+// wavefront restoration methods (according to functions' names)
+typedef struct{
+	WF_stat stat[WF_bench_size];
+	int Znum; // amount of Zernike coefficients
+	int measurements; // amount of measurements: std=sqrt(ZC_sum2/measurements-(ZC_sum/measurements)^2)
+} WF_benchmark;
+
+void free_wavefront(wavefront **F);
+double *calc_surface(int sz, double *Zidxs, int znum);
+wavefront *calc_surfaceR(int sz, double *Zidxs, int znum);
+wavefront *calc_surfaceRS(int sz, polar *points, double *Zidxs, int znum);
+polar *calc_BTA_Hpoints(int *sz);
+void free_bench(WF_benchmark **bench);
+void do_benchmark(double *Zidxs0, int znum0);
+#endif // __BENCHMARK_H__

cmdlnopts.c

@@ -0,0 +1,97 @@
+/*
+ * cmdlnopts.c - the only function that parse cmdln args and returns glob parameters
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+#include <assert.h>
+#include <stdio.h>
+#include <string.h>
+#include <strings.h>
+#include <math.h>
+#include "zernike.h"
+#include "cmdlnopts.h"
+#include "usefull_macros.h"
+
+#define RAD 57.2957795130823
+#define D2R(x) ((x) / RAD)
+#define R2D(x) ((x) * RAD)
+
+/*
+ * here are global parameters initialisation
+ */
+int help;
+glob_pars  G;
+
+int rewrite_ifexists = 0, // rewrite existing files == 0 or 1
+	verbose = 0; // each -v increments this value, e.g. -vvv sets it to 3
+//            DEFAULTS
+// default global parameters
+glob_pars const Gdefault = {
+	 .pixsize = 30.e-3
+	,.outfile = "wavefront"
+	,.imagesize = 1024
+	,.bench_maxP = 21  // Nmax = 5
+	,.bench_iter = 100
+};
+
+/*
+ * Define command line options by filling structure:
+ *	name	has_arg	flag	val		type		argptr			help
+*/
+myoption cmdlnopts[] = {
+	// set 1 to param despite of its repeating number:
+	{"help",	NO_ARGS,	NULL,	'h',	arg_int,	APTR(&help),		N_("show this help")},
+	{"verbose",	NO_ARGS,	NULL,	'v',	arg_int,	APTR(&verbose),		N_("be more and more verbose")},
+	{"rewrite",	NO_ARGS,	NULL,	'r',	arg_none,	APTR(&rewrite_ifexists),N_("rewrite existant files")},
+	{"gradients",NEED_ARG,	NULL,	'g',	arg_string,	APTR(&G.gradname),	N_("file with pre-computed gradients")},
+	{"pixsize",	NEED_ARG,	NULL,	'p',	arg_double,	APTR(&G.pixsize),	N_("CCD pixel size (if differs from FITS header)")},
+	{"zern-gen",NO_ARGS,	NULL,	'z',	arg_none,	APTR(&G.zerngen),	N_("generate FITS file with wavefront by given Zernike coefficients, in benchmark mode - use this coefficients instead of random")},
+	{"output",	NEED_ARG,	NULL,	'o',	arg_string,	APTR(&G.outfile),	N_("output file name (\"wavefront\" by default")},
+	{"benchmark",NO_ARGS,	NULL,	'b',	arg_none,	APTR(&G.benchmark),	N_("benchmark: test accuracy of different methods")},
+	{"imsize",	NEED_ARG,	NULL,	's',	arg_int,	APTR(&G.imagesize),	N_("size of output image (rectangular)")},
+	{"bench-maxP",NEED_ARG,	NULL,	'm',	arg_int,	APTR(&G.bench_maxP),N_("maximum amount of Zernike coefficients for benchmark")},
+	{"bench-iter",NEED_ARG,	NULL,	'i',	arg_int,	APTR(&G.bench_iter),N_("amount of iterations for benchmark (if -z is absent)")},
+	{"zprec",	NEED_ARG,	NULL,	'Z',	arg_double,	APTR(&Z_prec),		N_("minimal value of non-zero Zernike coefficients")},
+	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;
+	void *ptr;
+	ptr = memcpy(&G, &Gdefault, sizeof(G)); assert(ptr);
+	// format of help: "Usage: progname [args]\n"
+	change_helpstring(_("Usage: %s [args] <input files or Zernike coefficients>\n\n\tWhere args are:\n"));
+	// parse arguments
+	parseargs(&argc, &argv, cmdlnopts);
+	if(help) showhelp(-1, cmdlnopts);
+	if(argc > 0){
+		G.rest_pars_num = argc;
+		G.rest_pars = calloc(argc, sizeof(char*));
+		for (i = 0; i < argc; i++)
+			G.rest_pars[i] = strdup(argv[i]);
+	}
+	return &G;
+}
+

cmdlnopts.h

@@ -0,0 +1,49 @@
+/*
+ * cmdlnopts.h - comand line options for parceargs
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+#pragma once
+#ifndef __CMDLNOPTS_H__
+#define __CMDLNOPTS_H__
+
+#include "parseargs.h"
+
+
+typedef struct{
+	char **rest_pars;    // names of input files or zernike coefficients
+	int rest_pars_num;   // amount of files' names / coefficients
+	char *gradname;      // file with pre-computed gradients
+	double pixsize;      // CCD pixel size
+	int zerngen;         // generate FITS file with Zernike surface
+	char *outfile;       // output file name
+	int imagesize;       // output image size (1024 by default)
+	int benchmark;       // different methods benchmark
+	int bench_maxP;      // max ZC number for behchmark (default - 20)
+	int bench_iter;      // amount of iterations for benchmark (default - 100)
+} glob_pars;
+
+
+// default & global parameters
+extern glob_pars const Gdefault;
+extern int rewrite_ifexists, verbose;
+
+glob_pars *parse_args(int argc, char **argv);
+extern glob_pars  G;
+#endif // __CMDLNOPTS_H__

main.c

@@ -0,0 +1,201 @@
+/*
+ * Z-BTA_test.c - simple test for models of hartmannograms
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+#include <stdio.h>
+#include <math.h>
+#include <getopt.h>
+#include <stdarg.h>
+#include <string.h>
+#include "usefull_macros.h"
+#include "cmdlnopts.h"
+#include "zernike.h"
+#include "spots.h"
+#include "saveimg.h"
+#include "benchmark.h"
+
+
+void signals(int signo){
+	exit(signo);
+}
+
+int main(int argc, char **argv){
+	int i, j;
+	//double scale;
+	hartmann **images;
+	mirror *mir = NULL;
+	wavefront *comp;
+	double *Zidxs = NULL;
+	size_t L;
+	polar *coords = NULL;
+	point *grads = NULL;
+	initial_setup();
+	parse_args(argc, argv);
+	if(G.zerngen){ // user give his zernike coefficients
+		if(G.rest_pars_num < 1) ERRX(_("You should give at least one Zernike coefficient"));
+		Zidxs = MALLOC(double, G.rest_pars_num);
+		for(i = 0; i < G.rest_pars_num; ++i){
+			if(!str2double(&Zidxs[i], G.rest_pars[i])){
+				ERRX(_("Bad double number: %s!"), G.rest_pars[i]);
+			}
+		}
+	}
+	if(G.benchmark){
+		DBG("Benchmark");
+		do_benchmark(Zidxs, G.rest_pars_num);
+		FREE(Zidxs);
+		return 0;
+	}
+	if(G.zerngen){ // generate Zernike wavefront
+		DBG("Generate Zernike surface");
+		double *z = calc_surface(G.imagesize, Zidxs, G.rest_pars_num);
+		if(z){
+			writeimg(G.outfile, G.imagesize, mir, z);
+			FREE(z);
+		}
+		return 0;
+	}
+	if(!G.gradname){
+		DBG("Calculate wavefront by hartmannograms");
+		if(G.rest_pars_num < 2) ERR(_("You should give at least two file names: pre- and postfocal!"));
+		images = MALLOC(hartmann*, G.rest_pars_num + 1);
+		for(i = 0; i < G.rest_pars_num; ++i)
+			images[i] = read_spots(G.rest_pars[i]);
+		mir = calc_mir_coordinates(images);
+		getQ(mir);
+		calc_Hartmann_constant(mir);
+		spot_diagram *spot_dia = calc_spot_diagram(mir, mir->z07);
+		//printf("\nmirror's coordinates should be corrected to (%g, %g)\n",
+			//spot_dia->center.x, spot_dia->center.y);
+		printf("Projection of center on mirror shifted by (%g, %g), image shifted by (%g, %g)\n",
+			images[1]->center.x*HARTMANN_Z/distance, images[1]->center.y*HARTMANN_Z/distance,
+			images[1]->center.x*(FOCAL_R-HARTMANN_Z)/distance, images[1]->center.y*(FOCAL_R-HARTMANN_Z)/distance);
+		double tr = sqrt(images[1]->center.x*images[1]->center.x+images[1]->center.y*images[1]->center.y);
+		printf("Beam tilt is %g''\n", tr/distance*206265.);
+		calc_gradients(mir, spot_dia);
+		coords = MALLOC(polar, mir->spotsnum);
+		grads  = MALLOC(point, mir->spotsnum);
+		printf("\nGradients of aberrations (*1e-6):\nray#      x         y         dx         dy         R         phi\n");
+		for(j = 0, i = 0; i < 258; ++i){
+			if(!mir->got[i]) continue;
+			printf("%4d  %8.2f  %8.2f %10.6f %10.6f %10.2f %10.6f\n",
+				i, mir->spots[i].x, mir->spots[i].y,
+				mir->grads[i].x * 1e6, mir->grads[i].y * 1e6,
+				mir->pol_spots[i].r * MIR_R, mir->pol_spots[i].theta * 180. / M_PI);
+			memcpy(&coords[j], &mir->pol_spots[i], sizeof(polar));
+			memcpy(&grads[j], &mir->grads[i], sizeof(point));
+			grads[j].x *= 1e3;
+			grads[j].y *= 1e3;
+			j++;
+		}
+		L = mir->spotsnum;
+		//scale = mir->Rmax;
+		FREE(spot_dia);
+		FREE(mir);
+		for(i = 0; i < G.rest_pars_num; ++i)
+			h_free(&images[i]);
+		FREE(images);
+	}else{
+//		L = read_gradients(gradname, &coords, &grads, &scale);
+		return 0;
+	}
+/*
+	// spots information
+	printf(GREEN "\nSpots:\n" OLDCOLOR "\n      r\ttheta\tDx(mm/m)\tDy(mm/m)\n");
+	for(i = 0; i < L; i++){
+		printf("%8.1f%8.4f%8.4f%8.4f\n", coords[i].r, coords[i].theta,
+				grads[i].x*1000., grads[i].y*1000.);
+	}
+	*/
+	// gradients decomposition (Less squares, Zhao polinomials)
+	int Zsz, lastidx;
+	Zidxs = gradZdecomposeR(10, L, coords, grads, &Zsz, &lastidx);
+	//double *Zidxs = LS_gradZdecomposeR(10, L, coords, grads, &Zsz, &lastidx);
+//	Zidxs[1] = 0.;
+//	Zidxs[2] = 0.;
+	lastidx++;
+	printf("\nGradZ decompose, coefficients (%d):\n", lastidx);
+	for(i = 0; i < lastidx; i++) printf("%g, ", Zidxs[i]);
+	printf("\n\n");
+
+	int GridSize = 15;
+	comp = calc_surfaceR(GridSize, Zidxs, lastidx);
+	if(comp){
+		double zero_val = 0., N = 0., *zd = comp->zdata;
+		polar *rect = comp->coordinates;
+		for(j = GridSize-1; j > -1; j--){
+			int idx = j*GridSize;
+			for(i = 0; i < GridSize; i++,idx++){
+				if(zd[idx] > 1e-9){
+					zero_val += zd[idx];
+					N++;
+				}
+			}
+		}
+		zero_val /= N;
+		printf("zero: %g\nMirror surface:\n", zero_val);
+		for(j = GridSize-1; j > -1; j--){
+			int idx = j*GridSize;
+			for(i = 0; i < GridSize; i++,idx++){
+				if(rect[idx].r > 1. || rect[idx].r < 0.2){
+					printf("           ");
+				}else{
+					printf("%7.3f", zd[idx]  * MIR_R);
+				}
+			}
+			printf("\n");
+		}
+		free_wavefront(&comp);
+	}
+	// save matrix 100x100 with mirror deformations in Octave file format
+	FILE *f = fopen("file.out", "w");
+	if(f){
+		GridSize = 100;
+		comp = calc_surfaceR(GridSize, Zidxs, lastidx);
+		if(comp){
+			polar *rect = comp->coordinates;
+			double *zd = comp->zdata;
+			fprintf(f, "# generated by Z-BTA_test\n"
+				"# name: dev_matrix\n# type: matrix\n# rows: %d\n# columns: %d\n",
+					GridSize, GridSize);
+			for(j = 0; j < GridSize; j++){
+				int idx = j*GridSize;
+				for(i = 0; i < GridSize; i++,idx++){
+					if(rect[idx].r > 1. || rect[idx].r < 0.2)
+						fprintf(f, "0. ");
+					else
+						fprintf(f, "%g ", zd[idx]  * MIR_R);
+				}
+				fprintf(f,"\n");
+			}
+			free_wavefront(&comp);
+		}
+		fclose(f);
+	}
+	// save image
+	double *z = calc_surface(G.imagesize, Zidxs, lastidx);
+	if(z){
+		writeimg(G.outfile, G.imagesize, mir, z);
+		FREE(z);
+	}
+	FREE(Zidxs);
+	FREE(coords);
+	FREE(grads);
+	return 0;
+}

parseargs.c

@@ -0,0 +1,497 @@
+/*
+ * parseargs.c - parsing command line arguments & print help
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+#include <stdio.h>	// printf
+#include <getopt.h>	// getopt_long
+#include <stdlib.h>	// calloc, exit, strtoll
+#include <assert.h>	// assert
+#include <string.h> // strdup, strchr, strlen
+#include <strings.h>// strcasecmp
+#include <limits.h> // INT_MAX & so on
+#include <libintl.h>// gettext
+#include <ctype.h>	// isalpha
+#include "parseargs.h"
+#include "usefull_macros.h"
+
+char *helpstring = "%s\n";
+
+/**
+ * Change standard help header
+ * MAY consist ONE "%s" for progname
+ * @param str (i) - new format
+ */
+void change_helpstring(char *s){
+	int pcount = 0, scount = 0;
+	char *str = s;
+	// check `helpstring` and set it to default in case of error
+	for(; pcount < 2; str += 2){
+		if(!(str = strchr(str, '%'))) break;
+		if(str[1] != '%') pcount++; // increment '%' counter if it isn't "%%"
+		else{
+			str += 2; // pass next '%'
+			continue;
+		}
+		if(str[1] == 's') scount++; // increment "%s" counter
+	};
+	if(pcount > 1 || pcount != scount){ // amount of pcount and/or scount wrong
+		/// "îÅÐÒÁ×ÉÌØÎÙÊ ÆÏÒÍÁÔ ÓÔÒÏËÉ ÐÏÍÏÝÉ"
+		ERRX(_("Wrong helpstring!"));
+	}
+	helpstring = s;
+}
+
+/**
+ * Carefull atoll/atoi
+ * @param num (o) - returning value (or NULL if you wish only check number) - allocated by user
+ * @param str (i) - string with number must not be NULL
+ * @param t   (i) - T_INT for integer or T_LLONG for long long (if argtype would be wided, may add more)
+ * @return TRUE if conversion sone without errors, FALSE otherwise
+ */
+static bool myatoll(void *num, char *str, argtype t){
+	long long tmp, *llptr;
+	int *iptr;
+	char *endptr;
+	assert(str);
+	assert(num);
+	tmp = strtoll(str, &endptr, 0);
+	if(endptr == str || *str == '\0' || *endptr != '\0')
+		return FALSE;
+	switch(t){
+		case arg_longlong:
+			llptr = (long long*) num;
+			*llptr = tmp;
+		break;
+		case arg_int:
+		default:
+			if(tmp < INT_MIN || tmp > INT_MAX){
+				/// "ãÅÌÏÅ ×ÎÅ ÄÏÐÕÓÔÉÍÏÇÏ ÄÉÁÐÁÚÏÎÁ"
+				WARNX(_("Integer out of range"));
+				return FALSE;
+			}
+			iptr = (int*)num;
+			*iptr = (int)tmp;
+	}
+	return TRUE;
+}
+
+// the same as myatoll but for double
+// There's no NAN & INF checking here (what if they would be needed?)
+static bool myatod(void *num, const char *str, argtype t){
+	double tmp, *dptr;
+	float *fptr;
+	char *endptr;
+	assert(str);
+	tmp = strtod(str, &endptr);
+	if(endptr == str || *str == '\0' || *endptr != '\0')
+		return FALSE;
+	switch(t){
+		case arg_double:
+			dptr = (double *) num;
+			*dptr = tmp;
+		break;
+		case arg_float:
+		default:
+			fptr = (float *) num;
+			*fptr = (float)tmp;
+		break;
+	}
+	return TRUE;
+}
+
+/**
+ * Get index of current option in array options
+ * @param opt (i) - returning val of getopt_long
+ * @param options (i) - array of options
+ * @return index in array
+ */
+static int get_optind(int opt, myoption *options){
+	int oind;
+	myoption *opts = options;
+	assert(opts);
+	for(oind = 0; opts->name && opts->val != opt; oind++, opts++);
+	if(!opts->name || opts->val != opt) // no such parameter
+	showhelp(-1, options);
+	return oind;
+}
+
+/**
+ * reallocate new value in array of multiple repeating arguments
+ * @arg paptr - address of pointer to array (**void)
+ * @arg type - its type (for realloc)
+ * @return pointer to new (next) value
+ */
+void *get_aptr(void *paptr, argtype type){
+	int i = 1;
+	void **aptr = *((void***)paptr);
+	if(aptr){ // there's something in array
+		void **p = aptr;
+		while(*p++) ++i;
+	}
+	size_t sz = 0;
+	switch(type){
+		default:
+		case arg_none:
+			/// "îÅ ÍÏÇÕ ÉÓÐÏÌØÚÏ×ÁÔØ ÎÅÓËÏÌØËÏ ÐÁÒÁÍÅÔÒÏ× ÂÅÚ ÁÒÇÕÍÅÎÔÏ×!"
+			ERRX("Can't use multiple args with arg_none!");
+		break;
+		case arg_int:
+			sz = sizeof(int);
+		break;
+		case arg_longlong:
+			sz = sizeof(long long);
+		break;
+		case arg_double:
+			sz = sizeof(double);
+		break;
+		case arg_float:
+			sz = sizeof(float);
+		break;
+		case arg_string:
+			sz = 0;
+		break;
+	/*	case arg_function:
+			sz = sizeof(argfn *);
+		break;*/
+	}
+	aptr = realloc(aptr, (i + 1) * sizeof(void*));
+	*((void***)paptr) = aptr;
+	aptr[i] = NULL;
+	if(sz){
+		aptr[i - 1] = malloc(sz);
+	}else
+		aptr[i - 1] = &aptr[i - 1];
+	return aptr[i - 1];
+}
+
+
+/**
+ * Parse command line arguments
+ * ! If arg is string, then value will be strdup'ed!
+ *
+ * @param argc (io) - address of argc of main(), return value of argc stay after `getopt`
+ * @param argv (io) - address of argv of main(), return pointer to argv stay after `getopt`
+ * BE CAREFUL! if you wanna use full argc & argv, save their original values before
+ * 		calling this function
+ * @param options (i) - array of `myoption` for arguments parcing
+ *
+ * @exit: in case of error this function show help & make `exit(-1)`
+  */
+void parseargs(int *argc, char ***argv, myoption *options){
+	char *short_options, *soptr;
+	struct option *long_options, *loptr;
+	size_t optsize, i;
+	myoption *opts = options;
+	// check whether there is at least one options
+	assert(opts);
+	assert(opts[0].name);
+	// first we count how much values are in opts
+	for(optsize = 0; opts->name; optsize++, opts++);
+	// now we can allocate memory
+	short_options = calloc(optsize * 3 + 1, 1); // multiply by three for '::' in case of args in opts
+	long_options = calloc(optsize + 1, sizeof(struct option));
+	opts = options; loptr = long_options; soptr = short_options;
+	// in debug mode check the parameters are not repeated
+#ifdef EBUG
+	char **longlist = MALLOC(char*, optsize);
+	char *shortlist = MALLOC(char, optsize);
+#endif
+	// fill short/long parameters and make a simple checking
+	for(i = 0; i < optsize; i++, loptr++, opts++){
+		// check
+		assert(opts->name); // check name
+#ifdef EBUG
+		longlist[i] = strdup(opts->name);
+#endif
+		if(opts->has_arg){
+			assert(opts->type != arg_none); // check error with arg type
+			assert(opts->argptr);  // check pointer
+		}
+		if(opts->type != arg_none) // if there is a flag without arg, check its pointer
+			assert(opts->argptr);
+		// fill long_options
+		// don't do memcmp: what if there would be different alignment?
+		loptr->name		= opts->name;
+		loptr->has_arg	= (opts->has_arg < MULT_PAR) ? opts->has_arg : 1;
+		loptr->flag		= opts->flag;
+		loptr->val		= opts->val;
+		// fill short options if they are:
+		if(!opts->flag && opts->val){
+#ifdef EBUG
+			shortlist[i] = (char) opts->val;
+#endif
+			*soptr++ = opts->val;
+			if(loptr->has_arg) // add ':' if option has required argument
+				*soptr++ = ':';
+			if(loptr->has_arg == 2) // add '::' if option has optional argument
+				*soptr++ = ':';
+		}
+	}
+	// sort all lists & check for repeating
+#ifdef EBUG
+	int cmpstringp(const void *p1, const void *p2){
+		return strcmp(* (char * const *) p1, * (char * const *) p2);
+	}
+	int cmpcharp(const void *p1, const void *p2){
+		return (int)(*(char * const)p1 - *(char *const)p2);
+	}
+	qsort(longlist, optsize, sizeof(char *), cmpstringp);
+	qsort(shortlist,optsize, sizeof(char), cmpcharp);
+	char *prevl = longlist[0], prevshrt = shortlist[0];
+	for(i = 1; i < optsize; ++i){
+		if(longlist[i]){
+			if(prevl){
+				if(strcmp(prevl, longlist[i]) == 0) ERRX("double long arguments: --%s", prevl);
+			}
+			prevl = longlist[i];
+		}
+		if(shortlist[i]){
+			if(prevshrt){
+				if(prevshrt == shortlist[i]) ERRX("double short arguments: -%c", prevshrt);
+			}
+			prevshrt = shortlist[i];
+		}
+	}
+#endif
+	// now we have both long_options & short_options and can parse `getopt_long`
+	while(1){
+		int opt;
+		int oindex = 0, optind = 0; // oindex - number of option in argv, optind - number in options[]
+		if((opt = getopt_long(*argc, *argv, short_options, long_options, &oindex)) == -1) break;
+		if(opt == '?'){
+			opt = optopt;
+			optind = get_optind(opt, options);
+			if(options[optind].has_arg == NEED_ARG || options[optind].has_arg == MULT_PAR)
+				showhelp(optind, options); // need argument
+		}
+		else{
+			if(opt == 0 || oindex > 0) optind = oindex;
+			else optind = get_optind(opt, options);
+		}
+		opts = &options[optind];
+		if(opt == 0 && opts->has_arg == NO_ARGS) continue; // only long option changing integer flag
+		// now check option
+		if(opts->has_arg == NEED_ARG || opts->has_arg == MULT_PAR)
+			if(!optarg) showhelp(optind, options); // need argument
+		void *aptr;
+		if(opts->has_arg == MULT_PAR){
+			aptr = get_aptr(opts->argptr, opts->type);
+		}else
+			aptr = opts->argptr;
+		bool result = TRUE;
+		// even if there is no argument, but argptr != NULL, think that optarg = "1"
+		if(!optarg) optarg = "1";
+		switch(opts->type){
+			default:
+			case arg_none:
+				if(opts->argptr) *((int*)aptr) += 1; // increment value
+			break;
+			case arg_int:
+				result = myatoll(aptr, optarg, arg_int);
+			break;
+			case arg_longlong:
+				result = myatoll(aptr, optarg, arg_longlong);
+			break;
+			case arg_double:
+				result = myatod(aptr, optarg, arg_double);
+			break;
+			case arg_float:
+				result = myatod(aptr, optarg, arg_float);
+			break;
+			case arg_string:
+				result = (*((void**)aptr) = (void*)strdup(optarg));
+			break;
+			case arg_function:
+				result = ((argfn)aptr)(optarg);
+			break;
+		}
+		if(!result){
+			showhelp(optind, options);
+		}
+	}
+	*argc -= optind;
+	*argv += optind;
+}
+
+/**
+ * compare function for qsort
+ * first - sort by short options; second - sort arguments without sort opts (by long options)
+ */
+static int argsort(const void *a1, const void *a2){
+	const myoption *o1 = (myoption*)a1, *o2 = (myoption*)a2;
+	const char *l1 = o1->name, *l2 = o2->name;
+	int s1 = o1->val, s2 = o2->val;
+	int *f1 = o1->flag, *f2 = o2->flag;
+	// check if both options has short arg
+	if(f1 == NULL && f2 == NULL){ // both have short arg
+		return (s1 - s2);
+	}else if(f1 != NULL && f2 != NULL){ // both don't have short arg - sort by long
+		return strcmp(l1, l2);
+	}else{ // only one have short arg -- return it
+		if(f2) return -1; // a1 have short - it is 'lesser'
+		else return 1;
+	}
+}
+
+/**
+ * Show help information based on myoption->help values
+ * @param oindex (i)  - if non-negative, show only help by myoption[oindex].help
+ * @param options (i) - array of `myoption`
+ *
+ * @exit:  run `exit(-1)` !!!
+ */
+void showhelp(int oindex, myoption *options){
+	int max_opt_len = 0; // max len of options substring - for right indentation
+	const int bufsz = 255;
+	char buf[bufsz+1];
+	myoption *opts = options;
+	assert(opts);
+	assert(opts[0].name); // check whether there is at least one options
+	if(oindex > -1){ // print only one message
+		opts = &options[oindex];
+		printf("  ");
+		if(!opts->flag && isalpha(opts->val)) printf("-%c, ", opts->val);
+		printf("--%s", opts->name);
+		if(opts->has_arg == 1) printf("=arg");
+		else if(opts->has_arg == 2) printf("[=arg]");
+		printf("  %s\n", _(opts->help));
+		exit(-1);
+	}
+	// header, by default is just "progname\n"
+	printf("\n");
+	if(strstr(helpstring, "%s")) // print progname
+		printf(helpstring, __progname);
+	else // only text
+		printf("%s", helpstring);
+	printf("\n");
+	// count max_opt_len
+	do{
+		int L = strlen(opts->name);
+		if(max_opt_len < L) max_opt_len = L;
+	}while((++opts)->name);
+	max_opt_len += 14; // format: '-S , --long[=arg]' - get addition 13 symbols
+	opts = options;
+	// count amount of options
+	int N; for(N = 0; opts->name; ++N, ++opts);
+	if(N == 0) exit(-2);
+	// Now print all help (sorted)
+	opts = options;
+	qsort(opts, N, sizeof(myoption), argsort);
+	do{
+		int p = sprintf(buf, "  "); // a little indent
+		if(!opts->flag && opts->val) // .val is short argument
+			p += snprintf(buf+p, bufsz-p, "-%c, ", opts->val);
+		p += snprintf(buf+p, bufsz-p, "--%s", opts->name);
+		if(opts->has_arg == 1) // required argument
+			p += snprintf(buf+p, bufsz-p, "=arg");
+		else if(opts->has_arg == 2) // optional argument
+			p += snprintf(buf+p, bufsz-p, "[=arg]");
+		assert(p < max_opt_len); // there would be magic if p >= max_opt_len
+		printf("%-*s%s\n", max_opt_len+1, buf, _(opts->help)); // write options & at least 2 spaces after
+		++opts;
+	}while(--N);
+	printf("\n\n");
+	exit(-1);
+}
+
+/**
+ * get suboptions from parameter string
+ * @param str - parameter string
+ * @param opt - pointer to suboptions structure
+ * @return TRUE if all OK
+ */
+bool get_suboption(char *str, mysuboption *opt){
+	int findsubopt(char *par, mysuboption *so){
+		int idx = 0;
+		if(!par) return -1;
+		while(so[idx].name){
+			if(strcasecmp(par, so[idx].name) == 0) return idx;
+			++idx;
+		}
+		return -1; // badarg
+	}
+	bool opt_setarg(mysuboption *so, int idx, char *val){
+		mysuboption *soptr = &so[idx];
+		bool result = FALSE;
+		void *aptr = soptr->argptr;
+		switch(soptr->type){
+			default:
+			case arg_none:
+				if(soptr->argptr) *((int*)aptr) += 1; // increment value
+				result = TRUE;
+			break;
+			case arg_int:
+				result = myatoll(aptr, val, arg_int);
+			break;
+			case arg_longlong:
+				result = myatoll(aptr, val, arg_longlong);
+			break;
+			case arg_double:
+				result = myatod(aptr, val, arg_double);
+			break;
+			case arg_float:
+				result = myatod(aptr, val, arg_float);
+			break;
+			case arg_string:
+				result = (*((void**)aptr) = (void*)strdup(val));
+			break;
+			case arg_function:
+				result = ((argfn)aptr)(val);
+			break;
+		}
+		return result;
+	}
+	char *tok;
+	bool ret = FALSE;
+	char *tmpbuf;
+	tok = strtok_r(str, ":,", &tmpbuf);
+	do{
+		char *val = strchr(tok, '=');
+		int noarg = 0;
+		if(val == NULL){ // no args
+			val = "1";
+			noarg = 1;
+		}else{
+			*val++ = '\0';
+			if(!*val || *val == ':' || *val == ','){ // no argument - delimeter after =
+				val = "1"; noarg = 1;
+			}
+		}
+		int idx = findsubopt(tok, opt);
+		if(idx < 0){
+			/// "îÅÐÒÁ×ÉÌØÎÙÊ ÐÁÒÁÍÅÔÒ: %s"
+			WARNX(_("Wrong parameter: %s"), tok);
+			goto returning;
+		}
+		if(noarg && opt[idx].has_arg == NEED_ARG){
+			/// "%s: ÎÅÏÂÈÏÄÉÍ ÁÒÇÕÍÅÎÔ!"
+			WARNX(_("%s: argument needed!"), tok);
+			goto returning;
+		}
+		if(!opt_setarg(opt, idx, val)){
+			/// "îÅÐÒÁ×ÉÌØÎÙÊ ÁÒÇÕÍÅÎÔ \"%s\" ÐÁÒÁÍÅÔÒÁ \"%s\""
+			WARNX(_("Wrong argument \"%s\" of parameter \"%s\""), val, tok);
+			goto returning;
+		}
+	}while((tok = strtok_r(NULL, ":,", &tmpbuf)));
+	ret = TRUE;
+returning:
+	return ret;
+}

parseargs.h

@@ -0,0 +1,124 @@
+/*
+ * parseargs.h - headers for parsing command line arguments
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+#pragma once
+#ifndef __PARSEARGS_H__
+#define __PARSEARGS_H__
+
+#include <stdbool.h>// bool
+#include <stdlib.h>
+
+#ifndef TRUE
+	#define TRUE true
+#endif
+
+#ifndef FALSE
+	#define FALSE false
+#endif
+
+// macro for argptr
+#define APTR(x)   ((void*)x)
+
+// if argptr is a function:
+typedef  bool(*argfn)(void *arg);
+
+/*
+ * type of getopt's argument
+ * WARNING!
+ * My function change value of flags by pointer, so if you want to use another type
+ * make a latter conversion, example:
+ * 		char charg;
+ * 		int iarg;
+ * 		myoption opts[] = {
+ * 		{"value", 1, NULL, 'v', arg_int, &iarg, "char val"}, ..., end_option};
+ * 		..(parse args)..
+ * 		charg = (char) iarg;
+ */
+typedef enum {
+	arg_none = 0,	// no arg
+	arg_int,		// integer
+	arg_longlong,	// long long
+	arg_double,		// double
+	arg_float,		// float
+	arg_string,		// char *
+	arg_function	// parse_args will run function `bool (*fn)(char *optarg, int N)`
+} argtype;
+
+/*
+ * Structure for getopt_long & help
+ * BE CAREFUL: .argptr is pointer to data or pointer to function,
+ * 		conversion depends on .type
+ *
+ * ATTENTION: string `help` prints through macro PRNT(), bu default it is gettext,
+ * but you can redefine it before `#include "parseargs.h"`
+ *
+ * if arg is string, then value wil be strdup'ed like that:
+ * 		char *str;
+ * 		myoption opts[] = {{"string", 1, NULL, 's', arg_string, &str, "string val"}, ..., end_option};
+ *		*(opts[1].str) = strdup(optarg);
+ * in other cases argptr should be address of some variable (or pointer to allocated memory)
+ *
+ * NON-NULL argptr should be written inside macro APTR(argptr) or directly: (void*)argptr
+ *
+ * !!!LAST VALUE OF ARRAY SHOULD BE `end_option` or ZEROS !!!
+ *
+ */
+typedef enum{
+	NO_ARGS  = 0,  // first three are the same as in getopt_long
+	NEED_ARG = 1,
+	OPT_ARG  = 2,
+	MULT_PAR
+} hasarg;
+
+typedef struct{
+	// these are from struct option:
+	const char *name;		// long option's name
+	hasarg      has_arg;	// 0 - no args, 1 - nesessary arg, 2 - optionally arg, 4 - need arg & key can repeat (args are stored in null-terminated array)
+	int        *flag;		// NULL to return val, pointer to int - to set its value of val (function returns 0)
+	int         val;		// short opt name (if flag == NULL) or flag's value
+	// and these are mine:
+	argtype     type;		// type of argument
+	void       *argptr;		// pointer to variable to assign optarg value or function `bool (*fn)(char *optarg, int N)`
+	const char *help;		// help string which would be shown in function `showhelp` or NULL
+} myoption;
+
+/*
+ * Suboptions structure, almost the same like myoption
+ * used in parse_subopts()
+ */
+typedef struct{
+	const char *name;
+	hasarg      has_arg;
+	argtype     type;
+	void       *argptr;
+} mysuboption;
+
+// last string of array (all zeros)
+#define end_option {0,0,0,0,0,0,0}
+#define end_suboption {0,0,0,0}
+
+extern const char *__progname;
+
+void showhelp(int oindex, myoption *options);
+void parseargs(int *argc, char ***argv, myoption *options);
+void change_helpstring(char *s);
+bool get_suboption(char *str, mysuboption *opt);
+
+#endif // __PARSEARGS_H__

saveimg.c

@@ -0,0 +1,294 @@
+/*
+ * saveimg.c - functions to save data in png and FITS formats
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+#include "usefull_macros.h"
+#include "cmdlnopts.h" // for flag "-r", which will tell to rewrite existing file
+#include "saveimg.h"
+#include <png.h>
+#include <fitsio.h>
+#include <math.h>
+
+#define Stringify(x) #x
+#define OMP_FOR(x) _Pragma(Stringify(omp parallel for x))
+
+#ifndef THREAD_NUMBER
+	#define THREAD_NUMBER 4
+#endif
+
+#ifndef OMP_NUM_THREADS
+	#define OMP_NUM_THREADS THREAD_NUMBER
+#endif
+
+
+#define TRYFITS(f, ...)						\
+do{int status = 0; f(__VA_ARGS__, &status);	\
+	if (status){ ret = 0;					\
+		fits_report_error(stderr, status);	\
+		goto returning;}					\
+}while(0)
+#define WRITEKEY(...)						\
+do{ int status = 0;							\
+	fits_write_key(fp, __VA_ARGS__, &status);\
+	if(status)								\
+		fits_report_error(stderr, status);	\
+}while(0)
+
+/**
+ * Create filename as   outfile + number + "." + suffix
+ * 		number -- any number from 1 to 9999
+ * This function simply returns "outfile.suffix" when "-r" option is set
+ *
+ * @param outfile - file name
+ * @param suffix  - file suffix
+ * @return created filename or NULL
+ */
+char *createfilename(char* outfile, char* suffix){
+	FNAME();
+	struct stat filestat;
+	char buff[256], sfx[32];
+	if(suffix) snprintf(sfx, 31, ".%s", suffix);
+	else sfx[0] = 0; // no suffix
+	if(rewrite_ifexists){ // there was key "-f": simply return copy of filename
+		if(snprintf(buff, 255, "%s%s", outfile, sfx) < 1){
+			DBG("error");
+			return NULL;
+		}
+		DBG("(-r present) filename: %s", buff);
+		return strdup(buff);
+	}
+	int num;
+	if(!outfile) outfile = "";
+	for(num = 1; num < 10000; num++){
+		if(snprintf(buff, 255, "%s_%04d%s", outfile, num, sfx) < 1){
+			DBG("error");
+			return NULL;
+		}
+		if(stat(buff, &filestat)){ // || !S_ISREG(filestat.st_mode)) // OK, file not exists
+			DBG("filename: %s", buff);
+			return strdup(buff);
+		}
+	}
+	DBG("n: %s\n", buff);
+	WARN("Oops! All  numbers are busy or other error!");
+	return NULL;
+}
+
+typedef struct{
+	double *image;
+	double min;
+	double max;
+	double avr;
+	double std;
+} ImStat;
+
+static ImStat glob_stat;
+
+/**
+ * compute basics image statictics
+ * @param img - image data
+ * @param size - image size W*H
+ * @return
+ */
+void get_stat(double *img, size_t size){
+	FNAME();
+	if(glob_stat.image == img) return;
+	size_t i;
+	double pv, sum=0., sum2=0., sz=(double)size;
+	double max = -1., min = 1e15;
+	for(i = 0; i < size; i++){
+		pv = (double) *img++;
+		sum += pv;
+		sum2 += (pv * pv);
+		if(max < pv) max = pv;
+		if(min > pv) min = pv;
+	}
+	glob_stat.image = img;
+	glob_stat.avr = sum/sz;
+	glob_stat.std = sqrt(fabs(sum2/sz - glob_stat.avr*glob_stat.avr));
+	glob_stat.max = max;
+	glob_stat.min = min;
+	DBG("Image stat: max=%g, min=%g, avr=%g, std=%g", max, min, glob_stat.avr, glob_stat.std);
+}
+
+/**
+ * Save data to fits file
+ * @param filename - filename to save to
+ * @param sz  - image size: sz x sz
+ * @param imbox - image bounding box
+ * @data  image data
+ * @return 0 if failed
+ */
+int writefits(char *filename, size_t sz, mirror *mir, double *data){
+	FNAME();
+	long naxes[2] = {sz, sz};
+	static char* newname = NULL;
+	char buf[80];
+	int ret = 1;
+	time_t savetime = time(NULL);
+	fitsfile *fp;
+	if(!filename) return 0;
+	newname = realloc(newname, strlen(filename + 2));
+	sprintf(newname, "!%s", filename); // say cfitsio that file could be rewritten
+	TRYFITS(fits_create_file, &fp, newname);
+	TRYFITS(fits_create_img, fp, FLOAT_IMG, 2, naxes);
+	// FILE / Input file original name
+	WRITEKEY(TSTRING, "FILE", filename, "Input file original name");
+	WRITEKEY(TSTRING, "DETECTOR", "Wavefront model", "Detector model");
+	// IMAGETYP / object, flat, dark, bias, scan, eta, neon, push
+	WRITEKEY(TSTRING, "IMAGETYP", "object", "Image type");
+	// DATAMAX, DATAMIN / Max,min pixel value
+	WRITEKEY(TDOUBLE, "DATAMAX", &glob_stat.max, "Max data value");
+	WRITEKEY(TDOUBLE, "DATAMIN", &glob_stat.min, "Min data value");
+	// Some Statistics
+	WRITEKEY(TDOUBLE, "DATAAVR", &glob_stat.avr, "Average data value");
+	WRITEKEY(TDOUBLE, "DATASTD", &glob_stat.std, "Standart deviation of data value");
+	// DATE / Creation date (YYYY-MM-DDThh:mm:ss, UTC)
+	strftime(buf, 79, "%Y-%m-%dT%H:%M:%S", gmtime(&savetime));
+	WRITEKEY(TSTRING, "DATE", buf, "Creation date (YYYY-MM-DDThh:mm:ss, UTC)");
+	// DATE-OBS / DATE OF OBS.
+	WRITEKEY(TSTRING, "DATE-OBS", buf, "DATE OF OBS. (YYYY-MM-DDThh:mm:ss, local)");
+	// OBJECT  / Object name
+	WRITEKEY(TSTRING, "OBJECT", "Modeled object", "Object name");
+	// BINNING / Binning
+	WRITEKEY(TSTRING, "XBIN", "1", "Horizontal binning");
+	WRITEKEY(TSTRING, "YBIN", "1", "Vertical binning");
+	// PROG-ID / Observation program identifier
+	WRITEKEY(TSTRING, "PROG-ID", "Wavefront modeling", "Observation program identifier");
+	// AUTHOR / Author of the program
+	WRITEKEY(TSTRING, "AUTHOR", "Edward V. Emelianov", "Author of the program");
+	if(mir){
+		WRITEKEY(TDOUBLE, "ZBESTFOC", &mir->zbestfoc, "Z coordinate of best focus");
+		WRITEKEY(TDOUBLE, "Z70",      &mir->z07,      "Z coordinate of best PSF on 70\%");
+	}
+	TRYFITS(fits_write_img, fp, TDOUBLE, 1, sz * sz, data);
+	TRYFITS(fits_close_file, fp);
+returning:
+	return ret;
+}
+
+static uint8_t *rowptr = NULL;
+uint8_t *processRow(double *irow, size_t width, double min, double wd){
+	FREE(rowptr);
+	//float umax = ((float)(UINT16_MAX-1));
+	rowptr = MALLOC(uint8_t, width * 3);
+	OMP_FOR(shared(irow))
+	for(size_t i = 0; i < width; i++){
+		double gray = (irow[i] - min) / wd;
+		if(gray == 0.) continue;
+		int G = (int)(gray * 4.);
+		double x = 4.*gray - (double)G;
+		uint8_t *ptr = &rowptr[i*3];
+		uint8_t r = 0, g = 0, b = 0;
+		switch(G){
+			case 0:
+				g = (uint8_t)(255. * x + 0.5);
+				b = 255;
+			break;
+			case 1:
+				g = 255;
+				b = (uint8_t)(255. * (1. - x) + 0.5);
+			break;
+			case 2:
+				r = (uint8_t)(255. * x + 0.5);
+				g = 255;
+			break;
+			case 3:
+				r = 255;
+				g = (uint8_t)(255. * (1. - x) + 0.5);
+			break;
+			default:
+				r = 255;
+		}
+		ptr[0] = r; ptr[1] = g; ptr[2] = b;
+		//ptr[0] = ptr[1] = ptr[2] = gray*255;
+	}
+	return rowptr;
+}
+
+int writepng(char *filename, size_t sz, double *data){
+	FNAME();
+	int ret = 1;
+	FILE *fp = NULL;
+	png_structp pngptr = NULL;
+	png_infop infoptr = NULL;
+	double min = glob_stat.min, wd = glob_stat.max - min;
+	double *row;
+
+	if ((fp = fopen(filename, "w")) == NULL){
+		perror("Can't open png file");
+		ret = 0;
+		goto done;
+	}
+	if ((pngptr = png_create_write_struct(PNG_LIBPNG_VER_STRING,
+							NULL, NULL, NULL)) == NULL){
+		perror("Can't create png structure");
+		ret = 0;
+		goto done;
+	}
+	if ((infoptr = png_create_info_struct(pngptr)) == NULL){
+		perror("Can't create png info structure");
+		ret = 0;
+		goto done;
+	}
+	png_init_io(pngptr, fp);
+	png_set_compression_level(pngptr, 1);
+	png_set_IHDR(pngptr, infoptr, sz, sz, 8, PNG_COLOR_TYPE_RGB,//16, PNG_COLOR_TYPE_GRAY,
+				PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT,
+				PNG_FILTER_TYPE_DEFAULT);
+	png_write_info(pngptr, infoptr);
+	png_set_swap(pngptr);
+	size_t height = sz;
+	for(row = &data[sz*(height-1)]; height > 0; row -= sz, height--)
+		png_write_row(pngptr, (png_bytep)processRow(row, sz, min, wd));
+	png_write_end(pngptr, infoptr);
+done:
+	if(fp) fclose(fp);
+	if(pngptr) png_destroy_write_struct(&pngptr, &infoptr);
+	return ret;
+}
+
+/**
+ * Save data to image file[s] with format t
+ * @param name - filename prefix or NULL to save to "outXXXX.format"
+ * @param t - image[s] type[s]
+ * @param width, height - image size
+ * @param imbox - image bounding box (for FITS header)
+ * @param mirror - mirror parameters (for FITS header)
+ * @param data  image data
+ * @return number of saved images
+ */
+int writeimg(char *name, size_t sz, mirror *mir, double *data){
+	FNAME();
+	char *filename = NULL;
+	int ret = 0;
+	get_stat(data, sz*sz);
+	filename = createfilename(name, "fits");
+	if(filename){
+		ret = writefits(filename, sz, mir, data);
+		FREE(filename);
+	}
+	filename = createfilename(name, "png");
+	if(filename){
+		ret += writepng(filename, sz, data);
+		FREE(filename);
+	}
+	return ret;
+}

saveimg.h

@@ -0,0 +1,33 @@
+/*
+ * saveimg.h
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+#pragma once
+#ifndef __SAVEIMG_H__
+#define __SAVEIMG_H__
+
+#include <stdint.h>
+#include "spots.h"
+
+int writeimg(char *name, size_t sz, mirror *mir, double *data);
+
+char *createfilename(char* outfile, char* suffix);
+
+#endif // __SAVEIMG_H__

simple_list.c

@@ -0,0 +1,93 @@
+/*
+ * simple_list.c - simple one-direction list
+ *
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <err.h>
+
+#include "simple_list.h"
+
+#define MALLOC alloc_simple_list
+/**
+ * Memory allocation with control
+ * @param N - number of elements
+ * @param S - size of one element
+ * @return allocated memory or die
+ */
+static void *alloc_simple_list(size_t N, size_t S){
+    void *p = calloc(N, S);
+    if(!p) err(1, "malloc");
+    return p;
+}
+
+void (*listdata_free)(void *node_data) = NULL; // external function to free(listdata)
+
+#define FREE(arg) do{free(arg); arg = NULL;}while(0)
+
+/**
+ * add element v to list with root root (also this can be last element)
+ * @param root (io) - pointer to root (or last element) of list or NULL
+ * 						if *root == NULL, just created node will be placed there
+ * @param v         - data inserted
+ * @return pointer to created node
+ */
+List *list_add(List **root, void *v){
+	List *node, *last;
+	if((node = (List*) MALLOC(1, sizeof(List))) == 0)  return NULL; // allocation error
+	node->data = v; // insert data
+	if(root){
+		if(*root){ // there was root node - search last
+			last = *root;
+			while(last->next) last = last->next;
+			last->next = node; // insert pointer to new node into last element in list
+		}
+		if(!*root) *root = node;
+	}
+	return node;
+}
+
+/**
+ * set listdata_free()
+ * @param fn - function that freec memory of (node)
+ */
+void listfree_function(void (*fn)(void *node)){
+	listdata_free = fn;
+}
+
+/**
+ * remove all nodes in list
+ * @param root - pointer to root node
+ */
+void list_free(List **root){
+	List *node = *root, *next;
+	if(!root || !*root) return;
+	do{
+		next = node->next;
+		if(listdata_free)
+			listdata_free(node->data);
+		free(node);
+		node = next;
+	}while(node);
+	*root = NULL;
+}
+

simple_list.h

@@ -0,0 +1,39 @@
+/*
+ * simple_list.h - header file for simple list support
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+#pragma once
+#ifndef __SIMPLE_LIST_H__
+#define __SIMPLE_LIST_H__
+
+#define LIST_T(x) ((void*) x)
+
+typedef struct list_{
+	void *data;
+	struct list_ *next;
+} List;
+
+// add element v to list with root root (also this can be last element)
+List *list_add(List **root, void *v);
+// set listdata_free()
+void listfree_function(void (*fn)(void *node));
+void list_free(List **root);
+
+#endif // __SIMPLE_LIST_H__

spots.c

@@ -0,0 +1,556 @@
+/*
+ * spots.c
+ *
+ * Copyright 2015 Edward V. Emelianov <eddy@sao.ru, 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 2 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+#include <sys/stat.h> // stat
+#include <sys/mman.h> // mmap
+#include <stdio.h>
+#include <fcntl.h>
+#include <unistd.h>
+#include <err.h>
+#include <string.h>
+#include <assert.h>
+#include <math.h>
+#include "spots.h"
+#include "cmdlnopts.h"
+#include "usefull_macros.h"
+
+#define MM_TO_ARCSEC(x)  (x*206265./FOCAL_R)
+
+const double _4F = 4. * FOCAL_R;
+const double _2F = 2. * FOCAL_R;
+
+double distance = -1.; // distanse between pre- and postfocal images in millimeters
+
+void spots_free(List **spots){
+	listfree_function(free);
+	list_free(spots);
+	listfree_function(NULL);
+}
+
+/**
+ * Read spots-file, find center of hartmannogram & convert coordinates
+ * @param filename (i) - name of spots-file
+ * @return dynamically allocated hartmanogram structure
+ * COORDINATES ARE IN MILLIMETERS!!!
+ */
+hartmann *read_spots(char *filename){
+	assert(filename);
+	mmapbuf *M = NULL;
+	M = My_mmap(filename);
+	hartmann *H = MALLOC(hartmann, 1);
+	point *spots = H->spots;
+	uint8_t *got = H->got;
+	H->filename = strdup(filename);
+	char *pos = M->data, *epos = pos + M->len;
+	// readout list of spots
+	for(; pos && pos < epos; pos = strchr(pos+1, '\n')){
+		double x, y;
+		int id, a, b;
+		//if(3 != sscanf(pos, "%d %*s %*s %*s %*s %*s %*s %*s %*s %lf %lf", &id, &x, &y))
+		if(3 != sscanf(pos, "%d %*s %*s %*s %*s %lf %lf", &id, &x, &y))
+			continue;
+		a = id/100; b = id%100;
+		if(b <  32) id = a*32 + b; // main spots
+		else if(a == 2) id = 256;  // inner marker
+		else id = 257;             // outern marker
+		spots[id].x = x;
+#ifdef MIR_Y
+		spots[id].y = -y;
+#else
+		spots[id].y = y;
+#endif
+		got[id] = 1;
+	};
+	// get center: simply get center of each pair of opposite spots & calculate average
+	// IDs: xyy -- x(yy+16), yy=[0..15]
+	double xc = 0., yc = 0.;
+	int i, j, N = 0;
+	for(i = 0; i < 8; i++){ // circles
+		int id0 = i*32, id1 = id0 + 16;
+		for(j = 0; j < 16; j++, id0++, id1++){ // first half
+			if(!got[id0] || !got[id1]) continue; // there's no such pair
+			double c1, c2;
+			c1 = (spots[id0].x + spots[id1].x) / 2.;
+			c2 = (spots[id0].y + spots[id1].y) / 2.;
+			xc += c1;
+			yc += c2;
+			//printf("center: %g, %g\n", c1, c2);
+			++N;
+		}
+	}
+	xc /= (double) N;
+	yc /= (double) N;
+	//printf("Calculated center: %g, %g\n", xc, yc);
+	H->center.x = xc * G.pixsize;
+	H->center.y = yc * G.pixsize;
+	printf("get common center: (%.1f. %.1f)mm (pixsize=%g)\n", H->center.x, H->center.y, G.pixsize);
+	// convert coordinates to center & fill spots array of H
+	for(i = 0; i < 258; i++){
+		if(!got[i]) continue;
+		spots[i].x = (spots[i].x - xc) * G.pixsize;
+		spots[i].y = (spots[i].y - yc) * G.pixsize;
+		//printf("spot #%d: (%.1f, %.1f)mm\n", i, spots[i].x, spots[i].y);
+	}
+#ifdef MIR_Y
+	const double stp = M_PI/16., an0 = -M_PI_2 + stp/2., _2pi = 2.*M_PI;
+#else
+	const double stp = M_PI/16., an0 = +M_PI_2 - stp/2., _2pi = 2.*M_PI;
+#endif
+	double dmean = 0.;
+//	printf("RAYS' theta:\n");
+	for(i = 0; i < 32; i++){ // rays
+		int id = i, N=0;
+#ifdef MIR_Y
+		double sum = 0., refang = an0 + stp * (double)i;
+#else
+		double sum = 0., refang = an0 - stp * (double)i;
+#endif
+//		printf("ray %d: ", i);
+		for(j = 0; j < 8; j++, id+=32){ // circles
+			if(!got[id]) continue;
+			double theta = atan2(spots[id].y, spots[id].x);
+//			printf("%.5f,", theta);
+			sum += theta;
+			++N;
+		}
+		double meanang = sum/(double)N, delta = refang-meanang;
+		if(delta > _2pi) delta -= _2pi;
+		if(delta < 0.) delta += _2pi;
+		if(delta < 0.) delta += _2pi;
+//		printf("mean: %g, delta: %g\n", meanang, delta);
+		dmean += delta;
+	}
+	dmean /= 32.;
+	printf("MEAN delta: %g\n\n", dmean);
+	//if(!prefocal) dmean += M_PI;
+	double s,c;
+	sincos(dmean, &s, &c);
+//	printf("sin: %g, cos: %g\n",s,c);
+	// rotate data to ZERO
+	for(i = 0; i < 258; i++){
+		if(!got[i]) continue;
+		double x = spots[i].x, y = spots[i].y;
+		spots[i].x = x*c+y*s;
+		spots[i].y = -x*s+y*c;
+	}
+	My_munmap(M);
+	return H;
+}
+
+/**
+ * Calculate coordinates of points on mirror
+ * Modify coordinates of hartmannogram centers:
+ *  !!! the center beam on prefocal hartmannogram will have zero coordinates !!!
+ * @param H (io)  - array of pre- and postfocal H
+ * @param D      - distance from mirror top to prefocal image in millimeters
+ * @return mirror structure with coordinates & tan parameters
+ */
+mirror *calc_mir_coordinates(hartmann **H){
+	assert(H); assert(H[0]); assert(H[1]);
+	point *pre = H[0]->spots, *post = H[1]->spots, *prec = &H[0]->center, *postc = &H[1]->center;
+	mirror *mir = MALLOC(mirror, 1);
+	uint8_t *got_pre = H[0]->got, *got_post = H[1]->got, *mirgot = mir->got;
+	double *Z = MALLOC(double, 259); // 258 points + center beam
+	int i, iter;
+	point *tans = mir->tans, *spots = mir->spots, *mirc = &mir->center;
+	H[1]->center.x -= H[0]->center.x;
+	H[1]->center.y -= H[0]->center.y;
+	H[0]->center.x = 0.;
+	H[0]->center.y = 0.;
+	double dx = H[1]->center.x, dy = H[1]->center.y;
+	printf("H1 center: (%g, %g)\n", dx, dy);
+	double FCCnumerator = 0., FCCdenominator = 0.; // components of Fcc (Focus for minimal circle of confusion):
+	/*
+	 *         SUM (x_pre * tans_x + y_pre * tans_y)
+	 * Fcc = -----------------------------------------
+	 *             SUM (tanx_x^2 + tans_y^2)
+	 */
+	for(i = 0; i < 258; i++){ // check point pairs on pre/post
+		if(got_pre[i] && got_post[i]){
+			double tx, ty;
+			mirgot[i] = 1;
+			++mir->spotsnum;
+			//tx = (dx + post[i].x - pre[i].x)/distance;
+			tx = (post[i].x - pre[i].x)/distance;
+			tans[i].x = tx;
+			//ty = (dy + post[i].y - pre[i].y)/distance;
+			ty = (post[i].y - pre[i].y)/distance;
+			tans[i].y = ty;
+			FCCnumerator += pre[i].x*tx + pre[i].y*ty;
+			FCCdenominator += tx*tx + ty*ty;
+		}
+	}
+	mir->zbestfoc = -FCCnumerator/FCCdenominator;
+	double D = FOCAL_R - mir->zbestfoc;
+	printf("BEST (CC) focus: %g; D = %g\n", mir->zbestfoc, D);
+	printf("\nCalculate spots centers projections to mirror surface\n");
+	point tanc;
+	tanc.x = (dx + postc->x - prec->x)/distance;
+	tanc.y = (dy + postc->y - prec->y)/distance;
+	memcpy(&mir->tanc, &tanc, sizeof(point));
+	double Zerr = 10.;
+	/*
+	 * X = x_pre + (Z-D)*tans_x
+	 * Y = y_pre + (Z-D)*tans_y
+	 * Z = (X^2 + Y^2) / (4F)
+	 */
+	for(iter = 0; iter < 10 && Zerr > 1e-6; iter++){ // calculate points position on mirror
+		Zerr = 0.;
+		double x, y;
+		for(i = 0; i < 258; i++){
+			if(!mirgot[i]) continue;
+			x = pre[i].x + tans[i].x * (Z[i] - D);
+			y = pre[i].y + tans[i].y * (Z[i] - D);
+			spots[i].x = x;
+			spots[i].y = y;
+			double newZ = (x*x + y*y)/_4F;
+			double d = newZ - Z[i];
+			Zerr += d*d;
+			Z[i] = newZ;
+		}
+		x = prec->x + tanc.x * (Z[258] - D);
+		y = prec->y + tanc.y * (Z[258] - D);
+		mirc->x = x; mirc->y = y;
+		Z[258] = (x*x + y*y)/_4F;
+		printf("iteration %d, sum(dZ^2) = %g\n", iter, Zerr);
+	}
+	// now calculate polar coordinates relative to calculated center
+	double xc = mirc->x, yc = mirc->y;
+	polar *rtheta = mir->pol_spots;
+//	double Rmax = 0.;
+	for(i = 0; i < 258; i++){
+		if(!mirgot[i]) continue;
+		double x = spots[i].x - xc, y = spots[i].y - yc;
+		double r = sqrt(x*x + y*y);
+		rtheta[i].r = r;
+		rtheta[i].theta = atan2(y, x);
+//		if(Rmax < r) Rmax = r;
+	}
+//	Rmax = 3000.;
+	for(i = 0; i < 258; i++){
+		if(mirgot[i])
+	//		rtheta[i].r /= Rmax;
+		rtheta[i].r /= MIR_R;
+	}
+//	mir->Rmax = Rmax;
+	FREE(Z);
+	return mir;
+}
+
+/**
+ * Calculate Hartmann constant
+ * @param mir    (i) - filled mirror structure
+ * @param prefoc (i) - prefocal hartmannogram
+ * @return constant value
+ *
+ * Hartmann constant is
+ *      SUM(r_i^2*|F_i-F|)     200000
+ * T = -------------------- * --------
+ *          SUM(r_i)            F_m^2
+ *
+ *       where:
+ * F_m - mirror focus ratio
+ * r_i - radius of ith zone
+ * F_i - its focus (like in calc_mir_coordinates), counted from prefocal image
+ *       SUM(r_i * F_i)
+ * F = ------------------         - mean focus (counted from prefocal image)
+ *          SUM(r_i)
+ */
+double calc_Hartmann_constant(mirror *mir){
+	uint8_t *mirgot = mir->got;
+	int i, j;
+	point *tans = mir->tans, *pre = mir->prefocal;
+	polar *p_spots = mir->pol_spots;
+	double foc_i[8], r_i[8], Rsum = 0.;
+	double FCCnumerator = 0., FCCdenominator = 0.; // components of Fcc (Focus for minimal circle of confusion):
+	for(j = 0; j < 8; ++j){ // cycle by circles
+		double Rj = 0.;
+		int Nj = 0, idx = j*32;
+		for(i = 0; i < 32; ++i, ++idx){ // run through points on same circle
+			if(!mirgot[idx]) continue;
+			++Nj;
+			Rj += p_spots[idx].r;
+			double tx = tans[idx].x, ty = tans[idx].y;
+			FCCnumerator += pre[idx].x*tx + pre[idx].y*ty;
+			FCCdenominator += tx*tx + ty*ty;
+		}
+		foc_i[j] = -FCCnumerator/FCCdenominator;
+		r_i[j] = Rj / Nj;
+		Rsum += r_i[j];
+		printf("focus on R = %g is %g\n", MIR_R * r_i[j], foc_i[j]);
+	}
+	double F = 0.;
+	for(j = 0; j < 8; ++j){
+		F += foc_i[j] * r_i[j];
+	}
+	F /= Rsum;
+	double numerator = 0.;
+	for(j = 0; j < 8; ++j){
+		numerator += r_i[j]*r_i[j]*fabs(foc_i[j]-F);
+	}
+	printf("Mean focus is %g, numerator: %g, Rsum = %g\n", F, numerator, Rsum);
+	// multiply by MIR_R because r_i are normed by R
+	double T = MIR_R * 2e5/FOCAL_R/FOCAL_R*numerator/Rsum;
+	printf("\nHartmann value T = %g\n", T);
+	return T;
+}
+
+static int cmpdbl(const void * a, const void * b){
+	if (*(double*)a > *(double*)b)
+		return 1;
+	else return -1;
+}
+
+/**
+ * Calculate energy in circle of confusion
+ * @param mir    (i) - filled mirror structure
+ * @param hartm (i)  - zero-terminated list of hartmannograms
+ */
+void getQ(mirror *mir){
+	point *pre = mir->prefocal, *tans = mir->tans;
+	uint8_t *got = mir->got;
+	int i, j, N = mir->spotsnum, N03 = 0.3*N, N05 = 0.5*N, N07 = 0.7*N, N09 = 0.9*N;
+	double zstart = mir->zbestfoc - 3., z, zend = mir->zbestfoc + 3.01;
+	double *R = MALLOC(double, N);
+	double z03, z05, z07, z09, r03=1e3, r05=1e3, r07=1e3, r09=1e3;
+	printf("\nEnergy in circle of confusion\n");
+	printf("z        mean(R)'' std(R)''   R0.3''    R0.5''   R0.7''    R0.9''   Rmax''\n");
+	for(z = zstart; z < zend; z += 0.1){
+		j = 0;
+		double Rsum = 0., R2sum = 0.;
+		for(i = 0; i < 258; i++){
+			if(!got[i]) continue;
+			double x, y, R2, R1;
+			x = pre[i].x + tans[i].x * z;
+			y = pre[i].y + tans[i].y * z;
+			R2 = x*x + y*y;
+			R1 = sqrt(R2);
+			R[j] = R1;
+			R2sum += R2;
+			Rsum += R1;
+			++j;
+		}
+		qsort(R, N, sizeof(double), cmpdbl);
+		double R3 = R[N03], R5 = R[N05], R7 = R[N07], R9 = R[N09];
+/*		for(i = N-1; i; --i) if(R[i] < R7) break;
+		R7 = R[i];
+		for(i = N-1; i; --i) if(R[i] < R9) break;
+		R9 = R[i];*/
+		printf("%6.2f  %8.6f  %8.6f  %8.6f  %8.6f  %8.6f  %8.6f  %8.6f\n",
+			z, MM_TO_ARCSEC(Rsum/N), MM_TO_ARCSEC(sqrt((R2sum - Rsum*Rsum/N)/(N-1.))),
+			MM_TO_ARCSEC(R3), MM_TO_ARCSEC(R5),
+			MM_TO_ARCSEC(R7), MM_TO_ARCSEC(R9), MM_TO_ARCSEC(R[N-1]));
+		if(r03 > R3){
+			r03 = R3; z03 = z;
+		}
+		if(r05 > R5){
+			r05 = R5; z05 = z;
+		}
+		if(r07 > R7){
+			r07 = R7; z07 = z;
+		}
+		if(r09 > R9){
+			r09 = R9; z09 = z;
+		}
+	}
+	mir->z07 = z07;
+	printf("\ngot best values: z03=%g (r=%g''), z05=%g (r=%g''), z07=%g (r=%g''), z09=%g (r=%g'')\n",
+		z03, MM_TO_ARCSEC(r03), z05, MM_TO_ARCSEC(r05),
+		z07, MM_TO_ARCSEC(r07), z09, MM_TO_ARCSEC(r09));
+	printf("\nEnergy for z = %g\n   R,''      q(r)\n", z07);
+	for(j = 0, i = 0; i < 258; ++i){
+		if(!got[i]) continue;
+		double x, y;
+		x = pre[i].x + tans[i].x * z07;
+		y = pre[i].y + tans[i].y * z07;
+		R[j] = sqrt(x*x + y*y);
+		++j;
+	}
+	qsort(R, N, sizeof(double), cmpdbl);
+	for(i = 0; i < N; ++i){
+		printf("%8.6f  %8.6f\n", MM_TO_ARCSEC(R[i]), (1.+(double)i)/N);
+	}
+	FREE(R);
+}
+
+void h_free(hartmann **H){
+	if(!H || !*H) return;
+	FREE((*H)->filename);
+	FREE(*H);
+}
+
+/**
+ * Calculate spot diagram for given z value
+ * @param mir  (i) - filled mirror structure
+ * @param prefoc (i) - prefocal hartmannogram
+ * @return allocated structure of spot diagram
+ */
+spot_diagram *calc_spot_diagram(mirror *mir, double z){
+	int i;
+	spot_diagram *SD = MALLOC(spot_diagram, 1);
+	point *pre = mir->prefocal, *spots = SD->spots;
+	point *tans = mir->tans;
+	uint8_t *got = mir->got;
+	memcpy(SD->got, mir->got, 258);
+	SD->center.x = mir->prefoc_center.x + mir->tanc.x * z;
+	SD->center.y = mir->prefoc_center.y + mir->tanc.y * z;
+	printf("spots center: (%g, %g)\n", SD->center.x, SD->center.y);
+	printf("\nSpot diagram for z = %g (all in mm)\n ray#        x           y\n", z);
+	for(i = 0; i < 258; ++i){
+		if(!got[i]) continue;
+		spots[i].x = pre[i].x + tans[i].x * z;
+		spots[i].y = pre[i].y + tans[i].y * z;
+		printf("%4d   %10.6f   %10.6f\n", i, spots[i].x, spots[i].y);
+	}
+	return SD;
+}
+
+/**
+ * Calculate gradients of mirror surface aberrations
+ * @param mir       (io) - mirror (mir->grads will be filled by gradients)
+ * @param foc_spots  (i) - spot diagram for best focus
+ */
+void calc_gradients(mirror *mir, spot_diagram *foc_spots){
+	int i;
+	point *grads = mir->grads, *spots = foc_spots->spots;
+	printf("Common tilt: d/dx = %g mkm/m, d/dy = %g mkm/m\n",
+		-foc_spots->center.x / _2F * 1e6, -foc_spots->center.y / _2F * 1e6);
+	uint8_t *got = mir->got;
+	for(i = 0; i < 258; ++i){
+		if(!got[i]) continue;
+		grads[i].x = -(spots[i].x) / _2F;
+		grads[i].y = (spots[i].y) / _2F;
+	}
+}
+
+#if 0
+/**
+ *
+ * @param H      (i) - array of thwo hartmannograms (0 - prefocal, 1 - postfocal)
+ * @param coords (o) - array of gradients' coordinates on prefocal image (allocated here) - the same as H[0]->spots coordinates
+ * @param grads  (o) - gradients' array (allocated here)
+ * @param scale  (o) - scale of polar coordinate R (== Rmax)
+ * @return size of built arrays
+ */
+size_t get_gradients(hartmann *H[], polar **coords, point **grads, double *scale){
+	size_t Sz = 0, i, j, L0, L1;
+	assert(H); assert(H[0]); assert(H[1]);
+	spot *S0 = H[0]->spots, *S1;
+	List *CG_list = NULL, *curCG = NULL;
+	double Scale = G.pixsize * 1e-6 / distance / 2., S = 0.; // tg(2a)=dx/D -> a \approx dx/(2D)
+	L0 = H[0]->len;
+	L1 = H[1]->len;
+	/*
+	 * Both lists have sorted structure
+	 * but they can miss some points - that's why we should find exact points.
+	 * To store dinamycally data I use List
+	 */
+	for(i = 0; i < L0; ++i, ++S0){
+		int Id0 = S0->id;
+		for(S1 = H[1]->spots, j=0; j < L1; ++j, ++S1){
+			if(S1->id > Id0) break; // point with Id0 not found
+			if(S1->id == Id0){
+				CG *cg = MALLOC(CG, 1);
+				cg->id = Id0;
+				cg->x = S0->x; cg->y = S0->y;
+				cg->Dx = -(S1->x - S0->x)*Scale;
+				cg->Dy = -(S1->y - S0->y)*Scale;
+				Sz++;
+				if(CG_list)
+					curCG = list_add(&curCG, cg);
+				else
+					curCG = list_add(&CG_list, cg);
+				break;
+			}
+		}
+	}
+	polar *C = MALLOC(polar, Sz), *cptr = C;
+	point *G = MALLOC(point, Sz), *gptr = G;
+	curCG = CG_list;
+	for(i = 0; i < Sz; i++, cptr++, gptr++, curCG = curCG->next){
+		double x, y, length, R;
+		CG *cur = (CG*)curCG->data;
+		x = cur->x; y = cur->y;
+		R = sqrt(x*x + y*y);
+		cptr->r = R;
+		if(S < R) S = R; // find max R
+		cptr->theta = atan2(y, x);
+		x = cur->Dx; y = cur->Dy;
+		length = sqrt(1. + x*x + y*y); // length of vector for norm
+		gptr->x = x / length;
+		gptr->y = y / length;
+	}
+	cptr = C;
+	for(i = 0; i < Sz; i++, cptr++)
+		cptr->r /= S;
+	*scale = S;
+	*coords = C; *grads = G;
+	listfree_function(free);
+	list_free(&CG_list);
+	listfree_function(NULL);
+	return Sz;
+}
+
+/**
+ * Readout of gradients file calculated somewhere outside
+ *
+ * @param coords (o) - array with coordinates on mirror (in meters), ALLOCATED HERE!
+ * @param grads  (o) - array with gradients (meters per meters), ALLOCATED HERE!
+ * @param scale  (o) - scale on mirror (Rmax)
+ * @return - amount of points readed
+ */
+size_t read_gradients(char *gradname, polar **coords, point **grads, double *scale){
+	assert(gradname);
+	mmapbuf *M = NULL;
+	double Rmax = 0.;
+	M = My_mmap(gradname);
+	size_t L = 0;
+	char *pos = M->data, *epos = pos + M->len;
+	for(; pos && pos < epos; pos = strchr(pos+1, '\n')){
+		double x, y, gx, gy, R;
+		if(4 != sscanf(pos, "%lf %lf %lf %lf", &x, &y, &gx, &gy))
+			continue;
+		R = sqrt(x*x + y*y);
+		if(R > Rmax) Rmax = R;
+		L++;
+	}
+	printf("Found %zd points\n", L);
+	polar *C = MALLOC(polar, L), *cptr = C;
+	point *G = MALLOC(point, L), *gptr = G;
+	pos = M->data, epos = pos + M->len;
+	for(; pos && pos < epos; pos = strchr(pos+1, '\n')){
+		double x, y, gx, gy, R;
+		if(4 != sscanf(pos, "%lf %lf %lf %lf", &x, &y, &gx, &gy))
+			continue;
+		R = sqrt(x*x + y*y);
+		cptr->r = R / Rmax;
+		cptr->theta = atan2(y, x);
+		gptr->x = gx*1e6;
+		gptr->y = gy*1e6;
+		printf("%5.2f %5.2f %5.2f %5.2f\n",x, y, gx*1e6,gy*1e6);
+		cptr++, gptr++;
+	}
+	My_munmap(M);
+	*scale = Rmax;
+	*coords = C; *grads = G;
+	return L;
+}
+#endif

spots.h

@@ -0,0 +1,98 @@
+/*
+ * spots.h
+ *
+ * Copyright 2015 Edward V. Emelianov <eddy@sao.ru, 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 2 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+#pragma once
+#ifndef __SPOTS_H__
+#define __SPOTS_H__
+
+#include <stdint.h>
+#include "zernike.h"
+#include "simple_list.h"
+
+extern double pixsize;
+extern double distance;
+
+/*
+// spot center
+typedef struct{
+	int id;		// spot identificator
+	double x;	// coordinates of center
+	double y;
+} spot;
+*/
+
+// hartmannogram
+typedef struct{
+	char *filename;   // spots-filename
+	uint8_t got[258]; // == 1 if there's this spot on image
+	point center;     // coordinate of center
+	point spots[258]; // spotlist, markers have numbers 256 & 257
+	double focus;     // focus value for given image
+} hartmann;
+
+// mirror
+typedef struct{
+	uint8_t got[258]; // == 1 if there's this spot on both pre-focal & post-focal images
+	int spotsnum;     // total amount of spots used
+	point center;     // x&y coordinates of center beam
+	point tanc;       // tangent of center beam
+	point tans[258];  // tangents of beams == (r_post-r_pre)/d
+	point spots[258]; // x&y coordinates of spots on mirror surface (z=[x^2+y^2]/4f)
+	polar pol_spots[258]; // polar coordinates of spots on mirror surface according to center (norm by mirror R)
+	point grads[258]; // gradients to mirror surface
+	point prefocal[258];// coordinates of spots on prefocal image[s] (mean for many images)
+	point prefoc_center;// coordinate of prefocal hartmannogram[s] center
+	point postfocal[258]; // -//- on postfocal image[s]
+	double zbestfoc;  // Z-coordinate (from pre-focal image) of best focus for minimal circle of confusion
+	double z07;       // Z of best focus for q=0.7
+//	double Rmax;      // max R of spot (@ which pol_spots.r == 1.)
+} mirror;
+
+// spot diagram
+typedef struct{
+	uint8_t got[258];   // the same as above
+	point center;       // center beam's coordinate
+	point spots[258];   // spots
+	double z;           // z from prefocal image
+}spot_diagram;
+
+// gradients structure: point coordinates & gradient components
+typedef struct{
+	int id;
+	double x;
+	double y;
+	double Dx;
+	double Dy;
+} CG;
+
+hartmann *read_spots(char *filename);
+void h_free(hartmann **H);
+mirror *calc_mir_coordinates(hartmann **H);
+void getQ(mirror *mir);
+double calc_Hartmann_constant(mirror *mir);
+spot_diagram *calc_spot_diagram(mirror *mir, double z);
+void calc_gradients(mirror *mir, spot_diagram *foc_spots);
+
+/*
+size_t get_gradients(hartmann *H[], polar **coords, point **grads, double *scale);
+size_t read_gradients(char *gradname, polar **coords, point **grads, double *scale);
+*/
+#endif // __SPOTS_H__

usefull_macros.c

@@ -0,0 +1,360 @@
+/*
+ * usefull_macros.h - a set of usefull functions: memory, color etc
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+#include "usefull_macros.h"
+
+/**
+ * function for different purposes that need to know time intervals
+ * @return double value: time in seconds
+ */
+double dtime(){
+	double t;
+	struct timeval tv;
+	gettimeofday(&tv, NULL);
+	t = tv.tv_sec + ((double)tv.tv_usec)/1e6;
+	return t;
+}
+
+/******************************************************************************\
+ *                          Coloured terminal
+\******************************************************************************/
+int globErr = 0; // errno for WARN/ERR
+
+// pointers to coloured output printf
+int (*red)(const char *fmt, ...);
+int (*green)(const char *fmt, ...);
+int (*_WARN)(const char *fmt, ...);
+
+/*
+ * format red / green messages
+ * name: r_pr_, g_pr_
+ * @param fmt ... - printf-like format
+ * @return number of printed symbols
+ */
+int r_pr_(const char *fmt, ...){
+	va_list ar; int i;
+	printf(RED);
+	va_start(ar, fmt);
+	i = vprintf(fmt, ar);
+	va_end(ar);
+	printf(OLDCOLOR);
+	return i;
+}
+int g_pr_(const char *fmt, ...){
+	va_list ar; int i;
+	printf(GREEN);
+	va_start(ar, fmt);
+	i = vprintf(fmt, ar);
+	va_end(ar);
+	printf(OLDCOLOR);
+	return i;
+}
+/*
+ * print red error/warning messages (if output is a tty)
+ * @param fmt ... - printf-like format
+ * @return number of printed symbols
+ */
+int r_WARN(const char *fmt, ...){
+	va_list ar; int i = 1;
+	fprintf(stderr, RED);
+	va_start(ar, fmt);
+	if(globErr){
+		errno = globErr;
+		vwarn(fmt, ar);
+		errno = 0;
+		globErr = 0;
+	}else
+		i = vfprintf(stderr, fmt, ar);
+	va_end(ar);
+	i++;
+	fprintf(stderr, OLDCOLOR "\n");
+	return i;
+}
+
+static const char stars[] = "****************************************";
+/*
+ * notty variants of coloured printf
+ * name: s_WARN, r_pr_notty
+ * @param fmt ... - printf-like format
+ * @return number of printed symbols
+ */
+int s_WARN(const char *fmt, ...){
+	va_list ar; int i;
+	i = fprintf(stderr, "\n%s\n", stars);
+	va_start(ar, fmt);
+	if(globErr){
+		errno = globErr;
+		vwarn(fmt, ar);
+		errno = 0;
+		globErr = 0;
+	}else
+		i = +vfprintf(stderr, fmt, ar);
+	va_end(ar);
+	i += fprintf(stderr, "\n%s\n", stars);
+	i += fprintf(stderr, "\n");
+	return i;
+}
+int r_pr_notty(const char *fmt, ...){
+	va_list ar; int i;
+	i = printf("\n%s\n", stars);
+	va_start(ar, fmt);
+	i += vprintf(fmt, ar);
+	va_end(ar);
+	i += printf("\n%s\n", stars);
+	return i;
+}
+
+/**
+ * Run this function in the beginning of main() to setup locale & coloured output
+ */
+void initial_setup(){
+	// setup coloured output
+	if(isatty(STDOUT_FILENO)){ // make color output in tty
+		red = r_pr_; green = g_pr_;
+	}else{ // no colors in case of pipe
+		red = r_pr_notty; green = printf;
+	}
+	if(isatty(STDERR_FILENO)) _WARN = r_WARN;
+	else _WARN = s_WARN;
+	// Setup locale
+	setlocale(LC_ALL, "");
+	setlocale(LC_NUMERIC, "C");
+#if defined GETTEXT_PACKAGE && defined LOCALEDIR
+	bindtextdomain(GETTEXT_PACKAGE, LOCALEDIR);
+	textdomain(GETTEXT_PACKAGE);
+#endif
+}
+
+/******************************************************************************\
+ *                                  Memory
+\******************************************************************************/
+/*
+ * safe memory allocation for macro ALLOC
+ * @param N - number of elements to allocate
+ * @param S - size of single element (typically sizeof)
+ * @return pointer to allocated memory area
+ */
+void *my_alloc(size_t N, size_t S){
+	void *p = calloc(N, S);
+	if(!p) ERR("malloc");
+	//assert(p);
+	return p;
+}
+
+/**
+ * Mmap file to a memory area
+ *
+ * @param filename (i) - name of file to mmap
+ * @return stuct with mmap'ed file or die
+ */
+mmapbuf *My_mmap(char *filename){
+	int fd;
+	char *ptr;
+	size_t Mlen;
+	struct stat statbuf;
+	/// "îÅ ÚÁÄÁÎÏ ÉÍÑ ÆÁÊÌÁ!"
+	if(!filename) ERRX(_("No filename given!"));
+	if((fd = open(filename, O_RDONLY)) < 0)
+	/// "îÅ ÍÏÇÕ ÏÔËÒÙÔØ %s ÄÌÑ ÞÔÅÎÉÑ"
+		ERR(_("Can't open %s for reading"), filename);
+	if(fstat (fd, &statbuf) < 0)
+	/// "îÅ ÍÏÇÕ ×ÙÐÏÌÎÉÔØ stat %s"
+		ERR(_("Can't stat %s"), filename);
+	Mlen = statbuf.st_size;
+	if((ptr = mmap (0, Mlen, PROT_READ, MAP_PRIVATE, fd, 0)) == MAP_FAILED)
+	/// "ïÛÉÂËÁ mmap"
+		ERR(_("Mmap error for input"));
+	/// "îÅ ÍÏÇÕ ÚÁËÒÙÔØ mmap'ÎÕÔÙÊ ÆÁÊÌ"
+	if(close(fd)) ERR(_("Can't close mmap'ed file"));
+	mmapbuf *ret = MALLOC(mmapbuf, 1);
+	ret->data = ptr;
+	ret->len = Mlen;
+	return  ret;
+}
+
+void My_munmap(mmapbuf *b){
+	if(munmap(b->data, b->len)){
+	/// "îÅ ÍÏÇÕ munmap"
+		ERR(_("Can't munmap"));
+	}
+	FREE(b);
+}
+
+
+/******************************************************************************\
+ *                          Terminal in no-echo mode
+\******************************************************************************/
+static struct termios oldt, newt; // terminal flags
+static int console_changed = 0;
+// run on exit:
+void restore_console(){
+	if(console_changed)
+		tcsetattr(STDIN_FILENO, TCSANOW, &oldt); // return terminal to previous state
+	console_changed = 0;
+}
+
+// initial setup:
+void setup_con(){
+	if(console_changed) return;
+	tcgetattr(STDIN_FILENO, &oldt);
+	newt = oldt;
+	newt.c_lflag &= ~(ICANON | ECHO);
+	if(tcsetattr(STDIN_FILENO, TCSANOW, &newt) < 0){
+		/// "îÅ ÍÏÇÕ ÎÁÓÔÒÏÉÔØ ËÏÎÓÏÌØ"
+		WARN(_("Can't setup console"));
+		tcsetattr(STDIN_FILENO, TCSANOW, &oldt);
+		signals(0); //quit?
+	}
+	console_changed = 1;
+}
+
+/**
+ * Read character from console without echo
+ * @return char readed
+ */
+int read_console(){
+	int rb;
+	struct timeval tv;
+	int retval;
+	fd_set rfds;
+	FD_ZERO(&rfds);
+	FD_SET(STDIN_FILENO, &rfds);
+	tv.tv_sec = 0; tv.tv_usec = 10000;
+	retval = select(1, &rfds, NULL, NULL, &tv);
+	if(!retval) rb = 0;
+	else {
+		if(FD_ISSET(STDIN_FILENO, &rfds)) rb = getchar();
+		else rb = 0;
+	}
+	return rb;
+}
+
+/**
+ * getchar() without echo
+ * wait until at least one character pressed
+ * @return character readed
+ */
+int mygetchar(){ // getchar() without need of pressing ENTER
+	int ret;
+	do ret = read_console();
+	while(ret == 0);
+	return ret;
+}
+
+
+/******************************************************************************\
+ *                              TTY with select()
+\******************************************************************************/
+static struct termio oldtty, tty; // TTY flags
+static int comfd = -1; // TTY fd
+
+// run on exit:
+void restore_tty(){
+	if(comfd == -1) return;
+	ioctl(comfd, TCSANOW, &oldtty ); // return TTY to previous state
+	close(comfd);
+	comfd = -1;
+}
+
+#ifndef BAUD_RATE
+#define BAUD_RATE B9600
+#endif
+// init:
+void tty_init(char *comdev){
+	DBG("\nOpen port...\n");
+	if ((comfd = open(comdev,O_RDWR|O_NOCTTY|O_NONBLOCK)) < 0){
+		WARN("Can't use port %s\n",comdev);
+		ioctl(comfd, TCSANOW, &oldtty); // return TTY to previous state
+		close(comfd);
+		signals(0); // quit?
+	}
+	DBG(" OK\nGet current settings... ");
+	if(ioctl(comfd,TCGETA,&oldtty) < 0){  // Get settings
+		/// "îÅ ÍÏÇÕ ÐÏÌÕÞÉÔØ ÎÁÓÔÒÏÊËÉ"
+		WARN(_("Can't get settings"));
+		signals(0);
+	}
+	tty = oldtty;
+	tty.c_lflag     = 0; // ~(ICANON | ECHO | ECHOE | ISIG)
+	tty.c_oflag     = 0;
+	tty.c_cflag     = BAUD_RATE|CS8|CREAD|CLOCAL; // 9.6k, 8N1, RW, ignore line ctrl
+	tty.c_cc[VMIN]  = 0;  // non-canonical mode
+	tty.c_cc[VTIME] = 5;
+	if(ioctl(comfd,TCSETA,&tty) < 0){
+		/// "îÅ ÍÏÇÕ ÕÓÔÁÎÏ×ÉÔØ ÎÁÓÔÒÏÊËÉ"
+		WARN(_("Can't set settings"));
+		signals(0);
+	}
+	DBG(" OK\n");
+}
+/**
+ * Read data from TTY
+ * @param buff (o) - buffer for data read
+ * @param length   - buffer len
+ * @return amount of readed bytes
+ */
+size_t read_tty(uint8_t *buff, size_t length){
+	ssize_t L = 0;
+	fd_set rfds;
+	struct timeval tv;
+	int retval;
+	FD_ZERO(&rfds);
+	FD_SET(comfd, &rfds);
+	tv.tv_sec = 0; tv.tv_usec = 50000; // wait for 50ms
+	retval = select(comfd + 1, &rfds, NULL, NULL, &tv);
+	if (!retval) return 0;
+	if(FD_ISSET(comfd, &rfds)){
+		if((L = read(comfd, buff, length)) < 1) return 0;
+	}
+	return (size_t)L;
+}
+
+int write_tty(uint8_t *buff, size_t length){
+	ssize_t L = write(comfd, buff, length);
+	if((size_t)L != length){
+		/// "ïÛÉÂËÁ ÚÁÐÉÓÉ!"
+		WARN("Write error!");
+		return 1;
+	}
+	return 0;
+}
+
+
+/**
+ * Safely convert data from string to double
+ *
+ * @param num (o) - double number read from string
+ * @param str (i) - input string
+ * @return 1 if success, 0 if fails
+ */
+int str2double(double *num, const char *str){
+	double res;
+	char *endptr;
+	if(!str) return 0;
+	res = strtod(str, &endptr);
+	if(endptr == str || *str == '\0' || *endptr != '\0'){
+		/// "îÅÐÒÁ×ÉÌØÎÙÊ ÆÏÒÍÁÔ ÞÉÓÌÁ double!"
+		WARNX("Wrong double number format!");
+		return 0;
+	}
+	if(num) *num = res; // you may run it like myatod(NULL, str) to test wether str is double number
+	return 1;
+}

usefull_macros.h

@@ -0,0 +1,134 @@
+/*
+ * usefull_macros.h - a set of usefull macros: memory, color etc
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+#pragma once
+#ifndef __USEFULL_MACROS_H__
+#define __USEFULL_MACROS_H__
+
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <sys/mman.h>
+#include <unistd.h>
+#include <string.h>
+#include <stdio.h>
+#include <stdarg.h>
+#include <errno.h>
+#include <err.h>
+#include <locale.h>
+#if defined GETTEXT_PACKAGE && defined LOCALEDIR
+/*
+ * GETTEXT
+ */
+#include <libintl.h>
+#define _(String)				gettext(String)
+#define gettext_noop(String)	String
+#define N_(String)				gettext_noop(String)
+#else
+#define _(String)				(String)
+#define N_(String)				(String)
+#endif
+#include <stdlib.h>
+#include <termios.h>
+#include <termio.h>
+#include <sys/time.h>
+#include <sys/types.h>
+#include <stdint.h>
+
+
+// unused arguments with -Wall -Werror
+#define _U_    __attribute__((__unused__))
+
+/*
+ * Coloured messages output
+ */
+#define RED			"\033[1;31;40m"
+#define GREEN		"\033[1;32;40m"
+#define OLDCOLOR	"\033[0;0;0m"
+
+#ifndef FALSE
+#define FALSE (0)
+#endif
+
+#ifndef TRUE
+#define TRUE (1)
+#endif
+
+/*
+ * ERROR/WARNING messages
+ */
+extern int globErr;
+extern void signals(int sig);
+#define ERR(...) do{globErr=errno; _WARN(__VA_ARGS__); signals(9);}while(0)
+#define ERRX(...) do{globErr=0; _WARN(__VA_ARGS__); signals(9);}while(0)
+#define WARN(...) do{globErr=errno; _WARN(__VA_ARGS__);}while(0)
+#define WARNX(...) do{globErr=0; _WARN(__VA_ARGS__);}while(0)
+
+/*
+ * print function name, debug messages
+ * debug mode, -DEBUG
+ */
+#ifdef EBUG
+	#define FNAME() fprintf(stderr, "\n%s (%s, line %d)\n", __func__, __FILE__, __LINE__)
+	#define DBG(...) do{fprintf(stderr, "%s (%s, line %d): ", __func__, __FILE__, __LINE__); \
+					fprintf(stderr, __VA_ARGS__);			\
+					fprintf(stderr, "\n");} while(0)
+#else
+	#define FNAME()	 do{}while(0)
+	#define DBG(...) do{}while(0)
+#endif //EBUG
+
+/*
+ * Memory allocation
+ */
+#define ALLOC(type, var, size)  type * var = ((type *)my_alloc(size, sizeof(type)))
+#define MALLOC(type, size) ((type *)my_alloc(size, sizeof(type)))
+#define FREE(ptr)  do{if(ptr){free(ptr); ptr = NULL;}}while(0)
+
+double dtime();
+
+// functions for color output in tty & no-color in pipes
+extern int (*red)(const char *fmt, ...);
+extern int (*_WARN)(const char *fmt, ...);
+extern int (*green)(const char *fmt, ...);
+void * my_alloc(size_t N, size_t S);
+void initial_setup();
+
+// mmap file
+typedef struct{
+	char *data;
+	size_t len;
+} mmapbuf;
+mmapbuf *My_mmap(char *filename);
+void My_munmap(mmapbuf *b);
+
+void restore_console();
+void setup_con();
+int read_console();
+int mygetchar();
+
+void restore_tty();
+void tty_init(char *comdev);
+size_t read_tty(uint8_t *buff, size_t length);
+int write_tty(uint8_t *buff, size_t length);
+
+int str2double(double *num, const char *str);
+
+#endif // __USEFULL_MACROS_H__

zern_private.h

@@ -0,0 +1,46 @@
+/*
+ * zern_private.h - private variables for zernike.c & zernikeR.c
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+#pragma once
+#ifndef __ZERN_PRIVATE_H__
+#define __ZERN_PRIVATE_H__
+
+#include <math.h>
+#include <string.h>
+#include <err.h>
+#include <stdio.h>
+
+#ifndef iabs
+#define iabs(a)  (((a)<(0)) ? (-a) : (a))
+#endif
+
+extern double *FK;
+
+// zernike.c
+void build_factorial();
+void free_rpow(double ***Rpow, int n);
+void build_rpow(int W, int H, int n, double **Rad, double ***Rad_pow);
+double **build_rpowR(int n, int Sz, polar *P);
+
+// zernike_annular.c
+polar *conv_r(polar *r0, int Sz);
+
+#endif // __ZERN_PRIVATE_H__

zernike.c

@@ -0,0 +1,494 @@
+/*
+ * zernike.c - wavefront decomposition using Zernike polynomials
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+/*
+ * Literature:
+@ARTICLE{2007OExpr..1518014Z,
+   author = {{Zhao}, C. and {Burge}, J.~H.},
+    title = "{Orthonormal vector polynomials in a unit circle Part I: basis set derived from gradients of Zernike polynomials}",
+  journal = {Optics Express},
+     year = 2007,
+   volume = 15,
+    pages = {18014},
+      doi = {10.1364/OE.15.018014},
+   adsurl = {http://adsabs.harvard.edu/abs/2007OExpr..1518014Z},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+@ARTICLE{2008OExpr..16.6586Z,
+   author = {{Zhao}, C. and {Burge}, J.~H.},
+    title = "{Orthonormal vector polynomials in a unit circle, Part II : completing the basis set}",
+  journal = {Optics Express},
+     year = 2008,
+    month = apr,
+   volume = 16,
+    pages = {6586},
+      doi = {10.1364/OE.16.006586},
+   adsurl = {http://adsabs.harvard.edu/abs/2008OExpr..16.6586Z},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+ *
+ * !!!ATTENTION!!! Axe Y points upside-down!
+ */
+
+
+#include "zernike.h"
+#include "zern_private.h"
+#include "usefull_macros.h"
+
+double *FK = NULL;
+/**
+ * Build pre-computed array of factorials from 1 to 100
+ */
+void build_factorial(){
+	double F = 1.;
+	int i;
+	if(FK) return;
+	FK = MALLOC(double, 100);
+	FK[0] = 1.;
+	for(i = 1; i < 100; i++)
+		FK[i] = (F *= (double)i);
+}
+
+double Z_prec = 1e-6; // precision of Zernike coefficients
+
+
+/**
+ * Convert polynomial order in Noll notation into n/m
+ * @param p (i) - order of Zernike polynomial in Noll notation
+ * @param N (o) - order of polynomial
+ * @param M (o) - angular parameter
+ */
+void convert_Zidx(int p, int *N, int *M){
+	int n = (int) floor((-1.+sqrt(1.+8.*p)) / 2.);
+	*M = (int)(2.0*(p - n*(n+1.)/2. - 0.5*(double)n));
+	*N = n;
+}
+
+/**
+ * Free array of R powers with power n
+ * @param Rpow (i) - array to free
+ * @param n        - power of Zernike polinomial for that array (array size = n+1)
+ */
+void free_rpow(double ***Rpow, int n){
+	int i, N = n+1;
+	for(i = 0; i < N; i++) FREE((*Rpow)[i]);
+	FREE(*Rpow);
+}
+
+/**
+ * Build two arrays: with R and its powers (from 0 to n inclusive)
+ * for cartesian quoter I of matrix with size WxH
+ * @param W, H        - size of initial matrix
+ * @param n           - power of Zernike polinomial (array size = n+1)
+ * @param Rad (o)     - NULL or array with R in quater I
+ * @param Rad_pow (o) - NULL or array with powers of R
+ */
+void build_rpow(int W, int H, int n, double **Rad, double ***Rad_pow){
+	double Rnorm = fmax((W - 1.) / 2., (H - 1.) / 2.);
+	int i,j, k, N = n+1, w = (W+1)/2, h = (H+1)/2, S = w*h;
+	double **Rpow = MALLOC(double*, N); // powers of R
+	Rpow[0] = MALLOC(double, S);
+	for(j = 0; j < S; j++) Rpow[0][j] = 1.; // zero's power
+	double *R = MALLOC(double, S);
+	// first - fill array of R
+	double xstart = (W%2) ? 0. : 0.5, ystart = (H%2) ? 0. : 0.5;
+	for(j = 0; j < h; j++){
+		double *pt = &R[j*w], x, ww = w;
+		for(x = xstart; x < ww; x++, pt++){
+			double yy = (j + ystart)/Rnorm, xx = x/Rnorm;
+			*pt = sqrt(xx*xx+yy*yy);
+		}
+	}
+	for(i = 1; i < N; i++){ // Rpow - is quater I of cartesian coordinates ('cause other are fully simmetrical)
+		Rpow[i] = MALLOC(double, S);
+		double *rp = Rpow[i], *rpo = Rpow[i-1];
+		for(j = 0; j < h; j++){
+			int idx = j*w;
+			double *pt = &rp[idx], *pto = &rpo[idx], *r = &R[idx];
+			for(k = 0; k < w; k++, pt++, pto++, r++)
+				*pt = (*pto) * (*r); // R^{n+1} = R^n * R
+		}
+	}
+	if(Rad) *Rad = R;
+	else FREE(R);
+	if(Rad_pow) *Rad_pow = Rpow;
+	else free_rpow(&Rpow, n);
+}
+
+/**
+ * Calculate value of Zernike polynomial on rectangular matrix WxH pixels
+ * Center of matrix will be zero point
+ * Scale will be set by max(W/2,H/2)
+ * @param n - order of polynomial (max: 100!)
+ * @param m - angular parameter of polynomial
+ * @param W - width of output array
+ * @param H - height of output array
+ * @param norm (o) - (if !NULL) normalize factor
+ * @return array of Zernike polynomials on given matrix
+ */
+double *zernfun(int n, int m, int W, int H, double *norm){
+	double Z = 0., *Zarr = NULL;
+	bool erparm = false;
+	if(W < 2 || H < 2)
+		errx(1, "Sizes of matrix must be > 2!");
+	if(n > 100)
+		errx(1, "Order of Zernike polynomial must be <= 100!");
+	if(n < 0) erparm = true;
+	if(n < iabs(m)) erparm = true; // |m| must be <= n
+	int d = n - m;
+	if(d % 2) erparm = true; // n-m must differ by a prod of 2
+	if(erparm)
+		errx(1, "Wrong parameters of Zernike polynomial (%d, %d)", n, m);
+	if(!FK) build_factorial();
+	double Xc = (W - 1.) / 2., Yc = (H - 1.) / 2.; // coordinate of circle's middle
+	int i, j, k, m_abs = iabs(m), iup = (n-m_abs)/2, w = (W+1)/2;
+	double *R, **Rpow;
+	build_rpow(W, H, n, &R, &Rpow);
+	int SS = W * H;
+	double ZSum = 0.;
+	// now fill output matrix
+	Zarr = MALLOC(double, SS); // output matrix W*H pixels
+	for(j = 0; j < H; j++){
+		double *Zptr = &Zarr[j*W];
+		double Ryd = fabs(j - Yc);
+		int Ry = w * (int)Ryd; // Y coordinate on R matrix
+		for(i = 0; i < W; i++, Zptr++){
+			Z = 0.;
+			double Rxd = fabs(i - Xc);
+			int Ridx = Ry + (int)Rxd; // coordinate on R matrix
+			if(R[Ridx] > 1.) continue; // throw out points with R>1
+			// calculate R_n^m
+			for(k = 0; k <= iup; k++){ // Sum
+				double z = (1. - 2. * (k % 2)) * FK[n - k]        //       (-1)^k * (n-k)!
+					/(//----------------------------------- -----   -------------------------------
+						FK[k]*FK[(n+m_abs)/2-k]*FK[(n-m_abs)/2-k] // k!((n+|m|)/2-k)!((n-|m|)/2-k)!
+					);
+				Z += z * Rpow[n-2*k][Ridx];  // *R^{n-2k}
+			}
+			// normalize
+			double eps_m = (m) ? 1. : 2.;
+			Z *= sqrt(2.*(n+1.) / M_PI / eps_m);
+			double theta = atan2(j - Yc, i - Xc);
+			// multiply to angular function:
+			if(m){
+				if(m > 0)
+					Z *= cos(theta*(double)m_abs);
+				else
+					Z *= sin(theta*(double)m_abs);
+			}
+			*Zptr = Z;
+			ZSum += Z*Z;
+		}
+	}
+	if(norm) *norm = ZSum;
+	// free unneeded memory
+	FREE(R);
+	free_rpow(&Rpow, n);
+	return Zarr;
+}
+
+/**
+ * Zernike polynomials in Noll notation for square matrix
+ * @param p - index of polynomial
+ * @other params are like in zernfun
+ * @return zernfun
+ */
+double *zernfunN(int p, int W, int H, double *norm){
+	int n, m;
+	convert_Zidx(p, &n, &m);
+	return zernfun(n,m,W,H,norm);
+}
+
+/**
+ * Zernike decomposition of image 'image' WxH pixels
+ * @param Nmax (i) - maximum power of Zernike polinomial for decomposition
+ * @param W, H (i) - size of image
+ * @param image(i) - image itself
+ * @param Zsz  (o) - size of Z coefficients array
+ * @param lastIdx (o) - (if !NULL) last non-zero coefficient
+ * @return array of Zernike coefficients
+ */
+double *Zdecompose(int Nmax, int W, int H, double *image, int *Zsz, int *lastIdx){
+	int i, SS = W*H, pmax, maxIdx = 0;
+	double *Zidxs = NULL, *icopy = NULL;
+	pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation
+	Zidxs = MALLOC(double, pmax);
+	icopy = MALLOC(double, W*H);
+	memcpy(icopy, image, W*H*sizeof(double)); // make image copy to leave it unchanged
+	*Zsz = pmax;
+	for(i = 0; i < pmax; i++){ // now we fill array
+		double norm, *Zcoeff = zernfunN(i, W, H, &norm);
+		int j;
+		double *iptr = icopy, *zptr = Zcoeff, K = 0.;
+		for(j = 0; j < SS; j++, iptr++, zptr++)
+			K += (*zptr) * (*iptr) / norm; // multiply matrixes to get coefficient
+		Zidxs[i] = K;
+		if(fabs(K) < Z_prec){
+			Zidxs[i] = 0.;
+			continue; // there's no need to substract values that are less than our precision
+		}
+		maxIdx = i;
+		iptr = icopy; zptr = Zcoeff;
+		for(j = 0; j < SS; j++, iptr++, zptr++)
+			*iptr -= K * (*zptr); // subtract composed coefficient to reduce high orders values
+		FREE(Zcoeff);
+	}
+	if(lastIdx) *lastIdx = maxIdx;
+	FREE(icopy);
+	return Zidxs;
+}
+
+/**
+ * Zernike restoration of image WxH pixels by Zernike polynomials coefficients
+ * @param Zsz  (i) - number of elements in coefficients array
+ * @param Zidxs(i) - array with Zernike coefficients
+ * @param W, H (i) - size of image
+ * @return restored image
+ */
+double *Zcompose(int Zsz, double *Zidxs, int W, int H){
+	int i, SS = W*H;
+	double *image = MALLOC(double, SS);
+	for(i = 0; i < Zsz; i++){ // now we fill array
+		double K = Zidxs[i];
+		if(fabs(K) < Z_prec) continue;
+		double *Zcoeff = zernfunN(i, W, H, NULL);
+		int j;
+		double *iptr = image, *zptr = Zcoeff;
+		for(j = 0; j < SS; j++, iptr++, zptr++)
+			*iptr += K * (*zptr); // add next Zernike polynomial
+		FREE(Zcoeff);
+	}
+	return image;
+}
+
+
+/**
+ * Components of Zj gradient without constant factor
+ * 		all parameters are as in zernfun
+ * @return array of gradient's components
+ */
+point *gradZ(int n, int m, int W, int H, double *norm){
+	point *gZ = NULL;
+	bool erparm = false;
+	if(W < 2 || H < 2)
+		errx(1, "Sizes of matrix must be > 2!");
+	if(n > 100)
+		errx(1, "Order of gradient of Zernike polynomial must be <= 100!");
+	if(n < 1) erparm = true;
+	if(n < iabs(m)) erparm = true; // |m| must be <= n
+	int d = n - m;
+	if(d % 2) erparm = true; // n-m must differ by a prod of 2
+	if(erparm)
+		errx(1, "Wrong parameters of  gradient of Zernike polynomial (%d, %d)", n, m);
+	if(!FK) build_factorial();
+	double Xc = (W - 1.) / 2., Yc = (H - 1.) / 2.; // coordinate of circle's middle
+	int i, j, k, m_abs = iabs(m), iup = (n-m_abs)/2, isum = (n+m_abs)/2, w = (W+1)/2;
+	double *R, **Rpow;
+	build_rpow(W, H, n, &R, &Rpow);
+	int SS = W * H;
+	// now fill output matrix
+	gZ = MALLOC(point, SS);
+	double ZSum = 0.;
+	for(j = 0; j < H; j++){
+		point *Zptr = &gZ[j*W];
+		double Ryd = fabs(j - Yc);
+		int Ry = w * (int)Ryd; // Y coordinate on R matrix
+		for(i = 0; i < W; i++, Zptr++){
+			double Rxd = fabs(i - Xc);
+			int Ridx = Ry + (int)Rxd; // coordinate on R matrix
+			double Rcurr = R[Ridx];
+			if(Rcurr > 1. || fabs(Rcurr) < DBL_EPSILON) continue; // throw out points with R>1
+			double theta = atan2(j - Yc, i - Xc);
+			// components of grad Zj:
+
+			// 1. Theta_j
+			double Tj = 1., costm, sintm;
+			sincos(theta*(double)m_abs, &sintm, &costm);
+			if(m){
+				if(m > 0)
+					Tj = costm;
+				else
+					Tj = sintm;
+			}
+
+			// 2. dTheta_j/Dtheta
+			double dTj = 0.;
+			if(m){
+				if(m < 0)
+					dTj = m_abs * costm;
+				else
+					dTj = -m_abs * sintm;
+			}
+			// 3. R_j & dR_j/dr
+			double Rj = 0., dRj = 0.;
+			for(k = 0; k <= iup; k++){
+				double rj = (1. - 2. * (k % 2)) * FK[n - k]        //       (-1)^k * (n-k)!
+					/(//----------------------------------- -----   -------------------------------
+						FK[k]*FK[isum-k]*FK[iup-k] // k!((n+|m|)/2-k)!((n-|m|)/2-k)!
+					);
+//DBG("rj = %g (n=%d, k=%d)\n",rj, n, k);
+				int kidx = n - 2*k;
+				Rj += rj * Rpow[kidx][Ridx];  // *R^{n-2k}
+				if(kidx > 0)
+					dRj += rj * kidx * Rpow[kidx - 1][Ridx];
+			}
+			// normalisation for Zernike
+			double eps_m = (m) ? 1. : 2., sq =  sqrt(2.*(double)(n+1) / M_PI / eps_m);
+			Rj *= sq, dRj *= sq;
+			// 4. sin/cos
+			double sint, cost;
+			sincos(theta, &sint, &cost);
+
+			// projections of gradZj
+			double TdR = Tj * dRj, RdT = Rj * dTj / Rcurr;
+			Zptr->x = TdR * cost - RdT * sint;
+			Zptr->y = TdR * sint + RdT * cost;
+			// norm factor
+			ZSum += Zptr->x * Zptr->x + Zptr->y * Zptr->y;
+		}
+	}
+	if(norm) *norm = ZSum;
+	// free unneeded memory
+	FREE(R);
+	free_rpow(&Rpow, n);
+	return gZ;
+}
+
+/**
+ * Build components of vector polynomial Sj
+ * 		all parameters are as in zernfunN
+ * @return Sj(n,m) on image points
+ */
+point *zerngrad(int p, int W, int H, double *norm){
+	int n, m;
+	convert_Zidx(p, &n, &m);
+	if(n < 1) errx(1, "Order of gradient Z must be > 0!");
+	int m_abs = iabs(m);
+	int i,j;
+	double K = 1., K1 = 1.;///sqrt(2.*n*(n+1.));
+	point *Sj = MALLOC(point, W*H);
+	point *Zj =  gradZ(n, m, W, H, NULL);
+	double Zsum = 0.;
+	if(n == m_abs || n < 3){ // trivial case: n = |m| (in case of n=2,m=0 n'=0 -> grad(0,0)=0
+		for(j = 0; j < H; j++){
+			int idx = j*W;
+			point *Sptr = &Sj[idx], *Zptr = &Zj[idx];
+			for(i = 0; i < W; i++, Sptr++, Zptr++){
+				Sptr->x = K1*Zptr->x;
+				Sptr->y = K1*Zptr->y;
+				Zsum += Sptr->x * Sptr->x + Sptr->y * Sptr->y;
+			}
+		}
+	}else{
+		K = sqrt(((double)n+1.)/(n-1.));
+		//K1 /= sqrt(2.);
+		// n != |m|
+		// I run gradZ() twice! But another variant (to make array of Zj) can meet memory lack
+		point *Zj_= gradZ(n-2, m, W, H, NULL);
+		for(j = 0; j < H; j++){
+			int idx = j*W;
+			point *Sptr = &Sj[idx], *Zptr = &Zj[idx], *Z_ptr = &Zj_[idx];
+			for(i = 0; i < W; i++, Sptr++, Zptr++, Z_ptr++){
+				Sptr->x = K1*(Zptr->x - K * Z_ptr->x);
+				Sptr->y = K1*(Zptr->y - K * Z_ptr->y);
+				Zsum += Sptr->x * Sptr->x + Sptr->y * Sptr->y;
+			}
+		}
+		FREE(Zj_);
+	}
+	FREE(Zj);
+	if(norm) *norm = Zsum;
+	return Sj;
+}
+
+/**
+ * Decomposition of image with normals to wavefront by Zhao's polynomials
+ * 		all like Zdecompose
+ * @return array of coefficients
+ */
+double *gradZdecompose(int Nmax, int W, int H, point *image, int *Zsz, int *lastIdx){
+	int i, SS = W*H, pmax, maxIdx = 0;
+	double *Zidxs = NULL;
+	point *icopy = NULL;
+	pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation
+	Zidxs = MALLOC(double, pmax);
+	icopy = MALLOC(point, SS);
+	memcpy(icopy, image, SS*sizeof(point)); // make image copy to leave it unchanged
+	*Zsz = pmax;
+	for(i = 1; i < pmax; i++){ // now we fill array
+		double norm;
+		point *dZcoeff = zerngrad(i, W, H, &norm);
+		int j;
+		point *iptr = icopy, *zptr = dZcoeff;
+		double K = 0.;
+		for(j = 0; j < SS; j++, iptr++, zptr++)
+			K += (zptr->x*iptr->x + zptr->y*iptr->y) / norm; // multiply matrixes to get coefficient
+		if(fabs(K) < Z_prec)
+			continue; // there's no need to substract values that are less than our precision
+		Zidxs[i] = K;
+		maxIdx = i;
+		iptr = icopy; zptr = dZcoeff;
+		for(j = 0; j < SS; j++, iptr++, zptr++){
+			iptr->x -= K * zptr->x; // subtract composed coefficient to reduce high orders values
+			iptr->y -= K * zptr->y;
+		}
+		FREE(dZcoeff);
+	}
+	if(lastIdx) *lastIdx = maxIdx;
+	FREE(icopy);
+	return Zidxs;
+}
+
+/**
+ * Restoration of image with normals Zhao's polynomials coefficients
+ *		all like Zcompose
+ * @return restored image
+ */
+point *gradZcompose(int Zsz, double *Zidxs, int W, int H){
+	int i, SS = W*H;
+	point *image = MALLOC(point, SS);
+	for(i = 1; i < Zsz; i++){ // now we fill array
+		double K = Zidxs[i];
+		if(fabs(K) < Z_prec) continue;
+		point *Zcoeff = zerngrad(i, W, H, NULL);
+		int j;
+		point *iptr = image, *zptr = Zcoeff;
+		for(j = 0; j < SS; j++, iptr++, zptr++){
+			iptr->x += K * zptr->x;
+			iptr->y += K * zptr->y;
+		}
+		FREE(Zcoeff);
+	}
+	return image;
+}
+
+double *convGradIdxs(double *gradIdxs, int Zsz){
+	double *idxNew = MALLOC(double, Zsz);
+	int i;
+	for(i = 1; i < Zsz; i++){
+		int n,m;
+		convert_Zidx(i, &n, &m);
+		int j = ((n+2)*(n+4) + m) / 2;
+		if(j >= Zsz) j = 0;
+		idxNew[i] = (gradIdxs[i] - sqrt((n+3.)/(n+1.))*gradIdxs[j]);
+	}
+	return idxNew;
+}

zernike.h

@@ -0,0 +1,92 @@
+/*
+ * zernike.h
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+#pragma once
+#ifndef __ZERNIKE_H__
+#define __ZERNIKE_H__
+
+#include <stdbool.h>
+#include <stdlib.h>
+
+// focal ratio of BTA in millimeters
+#define FOCAL_R     (24024.)
+// BTA primary mirror radius (in mm)
+#define MIR_R       (3025.)
+// Distance from mirror to hartmann mask
+#define HARTMANN_Z  (20017.)
+
+/*************** Data structures & typedefs ***************/
+// point coordinates
+typedef struct{
+	double x,y;
+} point;
+
+typedef struct{
+	double r,theta;
+} polar;
+
+// 2D array
+typedef struct{
+	double **data;
+	size_t len; // size of 1D arrays
+	size_t num; // number of 1D arrays
+}_2D;
+
+extern double Z_prec; // precision of Zernike coefficients
+
+#ifndef DBL_EPSILON
+#define DBL_EPSILON        2.2204460492503131e-16
+#endif
+
+/*************** Base functions ***************/
+void convert_Zidx(int p, int *N, int *M);
+extern double Z_prec;
+
+/*************** Zernike on rectangular equidistant coordinate matrix ***************/
+double *zernfun(int n, int m, int W, int H, double *norm);
+double *zernfunN(int p, int W, int H, double *norm);
+
+double *Zdecompose(int Nmax, int W, int H, double *image, int *Zsz, int *lastIdx);
+double *Zcompose(int Zsz, double *Zidxs, int W, int H);
+
+double *gradZdecompose(int Nmax, int W, int H, point *image, int *Zsz, int *lastIdx);
+point *gradZcompose(int Zsz, double *Zidxs, int W, int H);
+double *convGradIdxs(double *gradIdxs, int Zsz);
+
+/*************** Zernike on a points set ***************/
+double *zernfunR(int n, int m, int Sz, polar *P, double *norm);
+double *zernfunNR(int p, int Sz, polar *P, double *norm);
+double *ZdecomposeR(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx);
+double *ZcomposeR(int Zsz, double *Zidxs, int Sz, polar *P);
+double *LS_decompose(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx);
+double *QR_decompose(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx);
+double *gradZdecomposeR(int Nmax, int Sz, polar *P,  point *grads, int *Zsz, int *lastIdx);
+double *LS_gradZdecomposeR(int Nmax, int Sz, polar *P,  point *grads, int *Zsz, int *lastIdx);
+point *gradZcomposeR(int Zsz, double *Zidxs, int Sz, polar *P);
+point *directGradZcomposeR(int Zsz, double *Zidxs, int Sz, polar *P);
+double *directGradZdecomposeR(int Nmax, int Sz, polar *P,  point *grads, int *Zsz, int *lastIdx);
+
+/*************** Annular Zernike ***************/
+_2D *ann_Z(int pmax, int Sz, polar *P, double **Norm);
+double *ann_Zcompose(int Zsz, double *Zidxs, int Sz, polar *P);
+double *ann_Zdecompose(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx);
+
+#endif // __ZERNIKE_H__

zernikeR.c

@@ -0,0 +1,687 @@
+/*
+ * zernikeR.c - Zernike decomposition for scattered points & annular aperture
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_blas.h>
+
+#include "usefull_macros.h"
+#include "zernike.h"
+#include "zern_private.h"
+
+/**
+ * Build array with R powers (from 0 to n inclusive)
+ * @param n     - power of Zernike polinomial (array size = n+1)
+ * @param Sz    - size of P array
+ * @param P (i) - polar coordinates of points
+ */
+double **build_rpowR(int n, int Sz, polar *P){
+	int i, j, N = n + 1;
+	double **Rpow = MALLOC(double*, N);
+	Rpow[0] = MALLOC(double, Sz);
+	for(i = 0; i < Sz; i++) Rpow[0][i] = 1.; // zero's power
+	for(i = 1; i < N; i++){ // Rpow - is quater I of cartesian coordinates ('cause other are fully simmetrical)
+		Rpow[i] = MALLOC(double, Sz);
+		double *rp = Rpow[i], *rpo = Rpow[i-1];
+		polar *p = P;
+		for(j = 0; j < Sz; j++, rp++, rpo++, p++){
+			*rp = (*rpo) * p->r;
+		}
+	}
+	return Rpow;
+}
+
+bool check_parameters(int n, int m, int Sz, polar *P){
+	bool erparm = false;
+	if(Sz < 3 || !P)
+		ERRX(_("Size of matrix must be > 2!"));
+	if(n > 100)
+		ERRX(_("Order of Zernike polynomial must be <= 100!"));
+	if(n < 0) erparm = true;
+	if(n < iabs(m)) erparm = true; // |m| must be <= n
+	if((n - m) % 2) erparm = true; // n-m must differ by a prod of 2
+	if(erparm)
+		ERRX(_("Wrong parameters of Zernike polynomial (%d, %d)"), n, m);
+	else
+		if(!FK) build_factorial();
+	return erparm;
+}
+
+/**
+ * Zernike function for scattering data
+ * @param n,m     - orders of polynomial
+ * @param Sz      - number of points
+ * @param P(i)    - array with points coordinates (polar, r<=1)
+ * @param norm(o) - (optional) norm coefficient
+ * @return dynamically allocated array with Z(n,m) for given array P
+ */
+double *zernfunR(int n, int m, int Sz, polar *P, double *norm){
+	if(check_parameters(n, m, Sz, P)) return NULL;
+	int j, k, m_abs = iabs(m), iup = (n-m_abs)/2;
+	double **Rpow = build_rpowR(n, Sz, P);
+	double ZSum = 0.;
+	// now fill output matrix
+	double *Zarr = MALLOC(double, Sz); // output matrix
+	double *Zptr = Zarr;
+	polar *p = P;
+	for(j = 0; j < Sz; j++, p++, Zptr++){
+		double Z = 0.;
+		if(p->r > 1.) continue; // throw out points with R>1
+		// calculate R_n^m
+		for(k = 0; k <= iup; k++){ // Sum
+			double z = (1. - 2. * (k % 2)) * FK[n - k]        //       (-1)^k * (n-k)!
+				/(//----------------------------------- -----   -------------------------------
+					FK[k]*FK[(n+m_abs)/2-k]*FK[(n-m_abs)/2-k] // k!((n+|m|)/2-k)!((n-|m|)/2-k)!
+				);
+			Z += z * Rpow[n-2*k][j];  // *R^{n-2k}
+		}
+		// normalize
+		double eps_m = (m) ? 1. : 2.;
+		Z *= sqrt(2.*(n+1.) / M_PI / eps_m  );
+		double m_theta = (double)m_abs * p->theta;
+		// multiply to angular function:
+		if(m){
+			if(m > 0)
+				Z *= cos(m_theta);
+			else
+				Z *= sin(m_theta);
+		}
+		*Zptr = Z;
+		ZSum += Z*Z;
+	}
+	if(norm) *norm = ZSum;
+	// free unneeded memory
+	free_rpow(&Rpow, n);
+	return Zarr;
+}
+
+/**
+ * Zernike polynomials in Noll notation for square matrix
+ * @param p - index of polynomial
+ * @other params are like in zernfunR
+ * @return zernfunR
+ */
+double *zernfunNR(int p, int Sz, polar *P, double *norm){
+	int n, m;
+	convert_Zidx(p, &n, &m);
+	return zernfunR(n,m,Sz,P,norm);
+}
+
+/**
+ * Zernike decomposition of WF with coordinates P
+ * @param Nmax (i)   - maximum power of Zernike polinomial for decomposition
+ * @param Sz, P(i)   - size (Sz) of points array (P)
+ * @param heights(i) - wavefront walues in points P
+ * @param Zsz  (o)   - size of Z coefficients array
+ * @param lastIdx(o) - (if !NULL) last non-zero coefficient
+ * @return array of Zernike coefficients
+ */
+double *ZdecomposeR(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx){
+	int i, pmax, maxIdx = 0;
+	double *Zidxs = NULL, *icopy = NULL;
+	pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation
+	Zidxs = MALLOC(double, pmax);
+	icopy = MALLOC(double, Sz);
+	memcpy(icopy, heights, Sz*sizeof(double)); // make image copy to leave it unchanged
+	*Zsz = pmax;
+	for(i = 0; i < pmax; i++){ // now we fill array
+		double norm, *Zcoeff = zernfunNR(i, Sz, P, &norm);
+		int j;
+		double *iptr = icopy, *zptr = Zcoeff, K = 0.;
+		for(j = 0; j < Sz; j++, iptr++, zptr++)
+			K += (*zptr) * (*iptr) / norm; // multiply matrixes to get coefficient
+		if(fabs(K) < Z_prec)
+			continue; // there's no need to substract values that are less than our precision
+		Zidxs[i] = K;
+		maxIdx = i;
+// Without these 3 lines wavefront decomposition on scattered and/or small amount of points would be wrong  ===================>
+		iptr = icopy; zptr = Zcoeff;
+		for(j = 0; j < Sz; j++, iptr++, zptr++)
+			*iptr -= K * (*zptr); // subtract composed coefficient to reduce high orders values
+//<=================== */
+		FREE(Zcoeff);
+	}
+	if(lastIdx) *lastIdx = maxIdx;
+	FREE(icopy);
+	return Zidxs;
+}
+
+/**
+ * Restoration of image in points P by Zernike polynomials coefficients
+ * @param Zsz  (i) - number of actual elements in coefficients array
+ * @param Zidxs(i) - array with Zernike coefficients
+ * @param Sz, P(i) - number (Sz) of points (P)
+ * @return restored image
+ */
+double *ZcomposeR(int Zsz, double *Zidxs, int Sz, polar *P){
+	int i;
+	double *image = MALLOC(double, Sz);
+	for(i = 0; i < Zsz; i++){ // now we fill array
+		double K = Zidxs[i];
+		if(fabs(K) < Z_prec) continue;
+		double *Zcoeff = zernfunNR(i, Sz, P, NULL);
+		int j;
+		double *iptr = image, *zptr = Zcoeff;
+		for(j = 0; j < Sz; j++, iptr++, zptr++)
+			*iptr += K * (*zptr); // add next Zernike polynomial
+		FREE(Zcoeff);
+	}
+	return image;
+}
+
+/**
+ * Prints out GSL matrix
+ * @param M (i) - matrix to print
+ */
+void print_matrix(gsl_matrix *M){
+	int x,y, H = M->size1, W = M->size2;
+	size_t T = M->tda;
+	printf("\nMatrix %dx%d:\n", W, H);
+	for(y = 0; y < H; y++){
+		double *ptr = &(M->data[y*T]);
+		printf("str %6d", y);
+		for(x = 0; x < W; x++, ptr++)
+			printf("%6.1f ", *ptr);
+		printf("\n");
+	}
+	printf("\n\n");
+}
+
+/**
+To try less-squares fit I decide after reading of
+@ARTICLE{2010ApOpt..49.6489M,
+   author = {{Mahajan}, V.~N. and {Aftab}, M.},
+    title = "{Systematic comparison of the use of annular and Zernike circle polynomials for annular wavefronts}",
+  journal = {\ao},
+     year = 2010,
+    month = nov,
+   volume = 49,
+    pages = {6489},
+      doi = {10.1364/AO.49.006489},
+   adsurl = {http://adsabs.harvard.edu/abs/2010ApOpt..49.6489M},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+*/
+/*
+ * n'th column of array m is polynomial of n'th degree,
+ * m'th row is m'th data point
+ *
+ * We fill matrix with polynomials by known datapoints coordinates,
+ * after that by less-square fit we get coefficients of decomposition
+ */
+double *LS_decompose(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx){
+	int pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation
+	if(Zsz) *Zsz = pmax;
+/*
+	(from GSL)
+	typedef struct	{
+		size_t size1;		// rows (height)
+		size_t size2;		// columns (width)
+		size_t tda;			// data width (aligned) - data stride
+		double * data;		// pointer to block->data (matrix data itself)
+		gsl_block * block;	// block with matrix data (block->size is size)
+		int owner;			// ==1 if matrix owns 'block' (then block will be freed with matrix)
+	} gsl_matrix;
+*/
+	// Now allocate matrix: its Nth row is equation for Nth data point,
+	// Mth column is Z_M
+	gsl_matrix *M = gsl_matrix_calloc(Sz, pmax);
+	// fill matrix with coefficients
+	int x,y;
+	size_t T = M->tda;
+	for(x = 0; x < pmax; x++){
+		double norm, *Zcoeff = zernfunNR(x, Sz, P, &norm), *Zptr = Zcoeff;
+		double *ptr = &(M->data[x]);
+		for(y = 0; y < Sz; y++, ptr+=T, Zptr++) // fill xth polynomial
+			*ptr = (*Zptr);
+		FREE(Zcoeff);
+	}
+
+	gsl_vector *ans = gsl_vector_calloc(pmax);
+
+	gsl_vector_view b = gsl_vector_view_array(heights, Sz);
+	gsl_vector *tau = gsl_vector_calloc(pmax); // min(size(M))
+	gsl_linalg_QR_decomp(M, tau);
+
+	gsl_vector *residual = gsl_vector_calloc(Sz);
+	gsl_linalg_QR_lssolve(M, tau, &b.vector, ans, residual);
+	double *ptr = ans->data;
+	int maxIdx = 0;
+	double *Zidxs = MALLOC(double, pmax);
+	for(x = 0; x < pmax; x++, ptr++){
+		if(fabs(*ptr) < Z_prec) continue;
+		Zidxs[x] = *ptr;
+		maxIdx = x;
+	}
+
+	gsl_matrix_free(M);
+	gsl_vector_free(ans);
+	gsl_vector_free(tau);
+	gsl_vector_free(residual);
+
+	if(lastIdx) *lastIdx = maxIdx;
+	return Zidxs;
+}
+
+/**
+ * Pseudo-annular Zernike polynomials
+ * They are just a linear composition of Zernike, made by Gram-Schmidt ortogonalisation
+ * Real Zernike koefficients restored by reverce transformation matrix
+ *
+ * Suppose we have a wavefront data in set of scattered points ${(x,y)}$, we want to
+ * find Zernike coefficients $z$ such that product of Zernike polynomials, $Z$, and
+ * $z$ give us that wavefront data with very little error (depending on $Z$ degree).
+ *
+ * Of cource, $Z$ isn't orthonormal on our basis, so we need to create some intermediate
+ * polynomials, $U$, which will be linear dependent from $Z$ (and this dependency
+ * should be strict and reversable, otherwise we wouldn't be able to reconstruct $z$
+ * from $u$): $U = Zk$. So, we have: $W = Uu + \epsilon$ and $W = Zz + \epsilon$.
+ *
+ * $U$ is orthonormal, so $U^T = U^{-1}$ and (unlike to $Z$) this reverce matrix
+ * exists and is unique. This mean that $u = W^T U = U^T W$.
+ * Coefficients matrix, $k$ must be inversable, so $Uk^{-1} = Z$, this mean that
+ * $z = uk^{-1}$.
+ *
+ * Our main goal is to find this matrix, $k$.
+ *
+ * 1. QR-decomposition
+ *  R. Navarro and J. Arines. Complete Modal Representation with Discrete Zernike
+ *  Polynomials - Critical Sampling in Non Redundant Grids. INTECH Open Access
+ *  Publisher, 2011.
+ *
+ * In this case non-orthogonal matrix $Z$ decomposed to orthogonal matrix $Q$ and
+ * right-triangular matrix $R$. In our case $Q$ is $U$ itself and $R$ is $k^{-1}$.
+ * QR-decomposition gives us an easy way to compute coefficient's matrix, $k$.
+ * Polynomials in $Q$ are linear-dependent from $Z$, each $n^{th}$ polynomial in $Q$
+ * depends from first $n$ polynomials in $Z$. So, columns of $R$ are coefficients
+ * for making $U$ from $Z$; rows in $R$ are coefficients for making $z$ from $u$.
+ *
+ * 2. Cholesky decomposition
+ * In this case for any non-orthogonal matrix with real values we have:
+ * $A^{T}A = LL^{T}$, where $L$ is left-triangular matrix.
+ * Orthogonal basis made by equation $U = A(L^{-1})^T$. And, as $A = UL^T$, we
+ * can reconstruct coefficients matrix: $k = (L^{-1})^T$.
+ */
+
+double *QR_decompose(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx){
+	int pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation
+	if(Sz < pmax) ERRX(_("Number of points must be >= number of polynomials!"));
+	if(Zsz) *Zsz = pmax;
+	int k, x,y;
+
+	//make new polar
+	polar *Pn = conv_r(P, Sz);
+	// Now allocate matrix: its Nth row is equation for Nth data point,
+	// Mth column is Z_M
+	gsl_matrix *M = gsl_matrix_calloc(Sz, pmax);
+	// Q-matrix (its first pmax columns is our basis)
+	gsl_matrix *Q = gsl_matrix_calloc(Sz, Sz);
+	// R-matrix (coefficients)
+	gsl_matrix *R = gsl_matrix_calloc(Sz, pmax);
+	// fill matrix with coefficients
+	size_t T = M->tda;
+	for(x = 0; x < pmax; x++){
+		double norm, *Zcoeff = zernfunNR(x, Sz, Pn, &norm), *Zptr = Zcoeff;
+		double *ptr = &(M->data[x]);
+		for(y = 0; y < Sz; y++, ptr+=T, Zptr++) // fill xth polynomial
+			*ptr = (*Zptr) / norm;
+		FREE(Zcoeff);
+	}
+
+	gsl_vector *tau = gsl_vector_calloc(pmax); // min(size(M))
+	gsl_linalg_QR_decomp(M, tau);
+
+	gsl_linalg_QR_unpack(M, tau, Q, R);
+//print_matrix(R);
+	gsl_matrix_free(M);
+	gsl_vector_free(tau);
+
+	double *Zidxs = MALLOC(double, pmax);
+	T = Q->tda;
+	size_t Tr = R->tda;
+
+	for(k = 0; k < pmax; k++){ // cycle by powers
+		double sumD = 0.;
+		double *Qptr = &(Q->data[k]);
+		for(y = 0; y < Sz; y++, Qptr+=T){ // cycle by points
+			sumD += heights[y] * (*Qptr);
+		}
+		Zidxs[k] = sumD;
+	}
+	gsl_matrix_free(Q);
+
+	// now restore Zernike coefficients
+	double *Zidxs_corr = MALLOC(double, pmax);
+	int maxIdx = 0;
+	for(k = 0; k < pmax; k++){
+		double sum = 0.;
+		double *Rptr = &(R->data[k*(Tr+1)]), *Zptr = &(Zidxs[k]);
+		for(x = k; x < pmax; x++, Zptr++, Rptr++){
+			sum += (*Zptr) * (*Rptr);
+		}
+		if(fabs(sum) < Z_prec) continue;
+		Zidxs_corr[k] = sum;
+		maxIdx = k;
+	}
+	gsl_matrix_free(R);
+	FREE(Zidxs);
+	FREE(Pn);
+	if(lastIdx) *lastIdx = maxIdx;
+	return Zidxs_corr;
+}
+
+
+/**
+ * Components of Zj gradient without constant factor
+ * @param n,m      - orders of polynomial
+ * @param Sz       - number of points
+ * @param P (i)    - coordinates of reference points
+ * @param norm (o) - norm factor or NULL
+ * @return array of gradient's components
+ */
+point *gradZR(int n, int m, int Sz, polar *P, double *norm){
+	if(check_parameters(n, m, Sz, P)) return NULL;
+	point *gZ = NULL;
+	int j, k, m_abs = iabs(m), iup = (n-m_abs)/2, isum = (n+m_abs)/2;
+	double **Rpow = build_rpowR(n, Sz, P);
+	// now fill output matrix
+	gZ = MALLOC(point, Sz);
+	point *Zptr = gZ;
+	double ZSum = 0.;
+	polar *p = P;
+	for(j = 0; j < Sz; j++, p++, Zptr++){
+		if(p->r > 1.) continue; // throw out points with R>1
+		double theta = p->theta;
+		// components of grad Zj:
+
+		// 1. Theta_j; 2. dTheta_j/Dtheta
+		double Tj = 1., dTj = 0.;
+		if(m){
+			double costm, sintm;
+			sincos(theta*(double)m_abs, &sintm, &costm);
+			if(m > 0){
+				Tj = costm;
+				dTj = -m_abs * sintm;
+			}else{
+				Tj = sintm;
+				dTj = m_abs * costm;
+			}
+		}
+
+		// 3. R_j & dR_j/dr
+		double Rj = 0., dRj = 0.;
+		for(k = 0; k <= iup; k++){
+			double rj = (1. - 2. * (k % 2)) * FK[n - k]        //       (-1)^k * (n-k)!
+				/(//----------------------------------- -----   -------------------------------
+					FK[k]*FK[isum-k]*FK[iup-k]                 // k!((n+|m|)/2-k)!((n-|m|)/2-k)!
+				);
+//DBG("rj = %g (n=%d, k=%d)\n",rj, n, k);
+			int kidx = n - 2*k;
+			Rj += rj * Rpow[kidx][j];  // *R^{n-2k}
+			if(kidx > 0)
+				dRj += rj * kidx * Rpow[kidx - 1][j]; // *(n-2k)*R^{n-2k-1}
+		}
+		// normalisation for Zernike
+		double eps_m = (m) ? 1. : 2., sq =  sqrt(2.*(double)(n+1) / M_PI / eps_m);
+		Rj *= sq, dRj *= sq;
+		// 4. sin/cos
+		double sint, cost;
+		sincos(theta, &sint, &cost);
+
+		// projections of gradZj
+		double TdR = Tj * dRj, RdT = Rj * dTj / p->r;
+		Zptr->x = TdR * cost - RdT * sint;
+		Zptr->y = TdR * sint + RdT * cost;
+		// norm factor
+		ZSum += Zptr->x * Zptr->x + Zptr->y * Zptr->y;
+	}
+	if(norm) *norm = ZSum;
+	// free unneeded memory
+	free_rpow(&Rpow, n);
+	return gZ;
+}
+
+/**
+ * Build components of vector polynomial Sj
+ * @param p        - index of polynomial
+ * @param Sz       - number of points
+ * @param P (i)    - coordinates of reference points
+ * @param norm (o) - norm factor or NULL
+ * @return Sj(n,m) on image points
+ */
+point *zerngradR(int p, int Sz, polar *P, double *norm){
+	int n, m, i;
+	convert_Zidx(p, &n, &m);
+	if(n < 1) ERRX(_("Order of gradient Z must be > 0!"));
+	int m_abs = iabs(m);
+	double Zsum, K = 1.;
+	point *Sj =  gradZR(n, m, Sz, P, &Zsum);
+	if(n != m_abs && n > 2){ // avoid trivial case: n = |m| (in case of n=2,m=0 n'=0 -> grad(0,0)=0
+		K = sqrt(((double)n+1.)/(n-1.));
+		Zsum = 0.;
+		point *Zj= gradZR(n-2, m, Sz, P, NULL);
+		point *Sptr = Sj, *Zptr = Zj;
+		for(i = 0; i < Sz; i++, Sptr++, Zptr++){
+			Sptr->x -= K * Zptr->x;
+			Sptr->y -= K * Zptr->y;
+			Zsum += Sptr->x * Sptr->x + Sptr->y * Sptr->y;
+		}
+		FREE(Zj);
+	}
+	if(norm) *norm = Zsum;
+	return Sj;
+}
+
+/**
+ * Decomposition of image with normals to wavefront by Zhao's polynomials
+ * by least squares method through LU-decomposition
+ * @param Nmax (i)   - maximum power of Zernike polinomial for decomposition
+ * @param Sz, P(i)   - size (Sz) of points array (P)
+ * @param grads(i) - wavefront gradients values in points P
+ * @param Zsz  (o)   - size of Z coefficients array
+ * @param lastIdx(o) - (if !NULL) last non-zero coefficient
+ * @return array of coefficients
+ */
+double *LS_gradZdecomposeR(int Nmax, int Sz, polar *P,  point *grads, int *Zsz, int *lastIdx){
+	int pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation
+	if(Zsz) *Zsz = pmax;
+	// Now allocate matrix: its Nth row is equation for Nth data point,
+	// Mth column is Z_M
+	int Sz2 = Sz*2, x, y;
+	gsl_matrix *M = gsl_matrix_calloc(Sz2, pmax);
+	// fill matrix with coefficients
+	size_t T = M->tda, S = T * Sz; // T is matrix stride, S - index of second coefficient
+	for(x = 1; x < pmax; x++){
+		double norm;
+		int n, m;
+		convert_Zidx(x, &n, &m);
+		point *dZcoeff = gradZR(n, m, Sz, P, &norm), *dZptr = dZcoeff;
+		//point *dZcoeff = zerngradR(x, Sz, P, &norm), *dZptr = dZcoeff;
+		double *ptr = &(M->data[x]);
+		// X-component is top part of resulting matrix, Y is bottom part
+		for(y = 0; y < Sz; y++, ptr+=T, dZptr++){ // fill xth polynomial
+			*ptr = dZptr->x;
+			ptr[S] = dZptr->y;
+		}
+		FREE(dZcoeff);
+	}
+
+	gsl_vector *ans = gsl_vector_calloc(pmax);
+	gsl_vector *b = gsl_vector_calloc(Sz2);
+	double *ptr = b->data;
+	for(x = 0; x < Sz; x++, ptr++, grads++){
+		// fill components of vector b like components of matrix M
+		*ptr = grads->x;
+		ptr[Sz] = grads->y;
+	}
+
+	gsl_vector *tau = gsl_vector_calloc(pmax);
+	gsl_linalg_QR_decomp(M, tau);
+
+	gsl_vector *residual = gsl_vector_calloc(Sz2);
+	gsl_linalg_QR_lssolve(M, tau, b, ans, residual);
+	ptr = &ans->data[1];
+	int maxIdx = 0;
+	double *Zidxs = MALLOC(double, pmax);
+	for(x = 1; x < pmax; x++, ptr++){
+		if(fabs(*ptr) < Z_prec) continue;
+		Zidxs[x] = *ptr;
+		maxIdx = x;
+	}
+
+	gsl_matrix_free(M);
+	gsl_vector_free(ans);
+	gsl_vector_free(b);
+	gsl_vector_free(tau);
+	gsl_vector_free(residual);
+
+	if(lastIdx) *lastIdx = maxIdx;
+	return Zidxs;
+}
+
+/**
+ * Decomposition of image with normals to wavefront by Zhao's polynomials
+ * @param Nmax (i)   - maximum power of Zernike polinomial for decomposition
+ * @param Sz, P(i)   - size (Sz) of points array (P)
+ * @param grads(i) - wavefront gradients values in points P
+ * @param Zsz  (o)   - size of Z coefficients array
+ * @param lastIdx(o) - (if !NULL) last non-zero coefficient
+ * @return array of coefficients
+ */
+double *gradZdecomposeR(int Nmax, int Sz, polar *P,  point *grads, int *Zsz, int *lastIdx){
+	int i, pmax, maxIdx = 0;
+	pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation
+	double *Zidxs = MALLOC(double, pmax);
+	point *icopy = MALLOC(point, Sz);
+	memcpy(icopy, grads, Sz*sizeof(point)); // make image copy to leave it unchanged
+	*Zsz = pmax;
+	for(i = 1; i < pmax; i++){ // now we fill array
+		double norm;
+		point *dZcoeff = zerngradR(i, Sz, P, &norm);
+		int j;
+		point *iptr = icopy, *zptr = dZcoeff;
+		double K = 0.;
+		for(j = 0; j < Sz; j++, iptr++, zptr++)
+			K += zptr->x*iptr->x + zptr->y*iptr->y; // multiply matrixes to get coefficient
+		K /= norm;
+		if(fabs(K) < Z_prec)
+			continue; // there's no need to substract values that are less than our precision
+		Zidxs[i] = K;
+		maxIdx = i;
+		iptr = icopy; zptr = dZcoeff;
+		for(j = 0; j < Sz; j++, iptr++, zptr++){
+			iptr->x -= K * zptr->x; // subtract composed coefficient to reduce high orders values
+			iptr->y -= K * zptr->y;
+		}
+		FREE(dZcoeff);
+	}
+	if(lastIdx) *lastIdx = maxIdx;
+	FREE(icopy);
+	return Zidxs;
+}
+
+/**
+ * Restoration of wavefront normals by given Zhao's polynomials coefficients
+ *		all like Zcompose, but `Zidxs` are Zhao's (not direct Zernike)!
+ * @return restored image
+ */
+point *gradZcomposeR(int Zsz, double *Zidxs, int Sz, polar *P){
+	int i;
+	point *image = MALLOC(point, Sz);
+	for(i = 1; i < Zsz; i++){ // now we fill array
+		double K = Zidxs[i];
+		if(fabs(K) < Z_prec) continue;
+		point *Zcoeff = zerngradR(i, Sz, P, NULL);
+		int j;
+		point *iptr = image, *zptr = Zcoeff;
+		for(j = 0; j < Sz; j++, iptr++, zptr++){
+			iptr->x += K * zptr->x;
+			iptr->y += K * zptr->y;
+		}
+		FREE(Zcoeff);
+	}
+	return image;
+}
+
+/**
+ * Direct gradients of Zernike polynomials
+ */
+point *directGradZcomposeR(int Zsz, double *Zidxs, int Sz, polar *P){
+	point *image = MALLOC(point, Sz);
+	for(int p = 1; p < Zsz; ++p){
+		int n, m;
+		double K =  Zidxs[p];
+		convert_Zidx(p, &n, &m);
+		point *Zcoeff =  gradZR(n, m, Sz, P, NULL);
+		point *iptr = image, *zptr = Zcoeff;
+		for(int j = 0; j < Sz; ++j, ++iptr, ++zptr){
+			iptr->x += K * zptr->x;
+			iptr->y += K * zptr->y;
+		}
+		FREE(Zcoeff);
+	}
+	return image;
+}
+
+/**
+ * Direct decomposition of image with normals to wavefront by Zernike gradients
+ * @param Nmax (i)   - maximum power of Zernike polinomial for decomposition
+ * @param Sz, P(i)   - size (Sz) of points array (P)
+ * @param grads(i) - wavefront gradients values in points P
+ * @param Zsz  (o)   - size of Z coefficients array
+ * @param lastIdx(o) - (if !NULL) last non-zero coefficient
+ * @return array of coefficients
+ */
+double *directGradZdecomposeR(int Nmax, int Sz, polar *P,  point *grads, int *Zsz, int *lastIdx){
+	int i, pmax, maxIdx = 0;
+	pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation
+	double *Zidxs = MALLOC(double, pmax);
+	point *icopy = MALLOC(point, Sz);
+	memcpy(icopy, grads, Sz*sizeof(point)); // make image copy to leave it unchanged
+	*Zsz = pmax;
+	for(i = 1; i < pmax; i++){ // now we fill array
+		double norm;
+		int n, m;
+		convert_Zidx(i, &n, &m);
+		point *dZcoeff = gradZR(n, m, Sz, P, &norm);
+		int j;
+		point *iptr = icopy, *zptr = dZcoeff;
+		double K = 0.;
+		for(j = 0; j < Sz; j++, iptr++, zptr++)
+			K += zptr->x*iptr->x + zptr->y*iptr->y;
+		K /= norm;
+		if(fabs(K) < Z_prec)
+			continue; // there's no need to substract values that are less than our precision
+		Zidxs[i] = K;
+		maxIdx = i;
+
+		iptr = icopy; zptr = dZcoeff;
+		for(j = 0; j < Sz; j++, iptr++, zptr++){
+			iptr->x -= K * zptr->x; // subtract composed coefficient to reduce high orders values
+			iptr->y -= K * zptr->y;
+		}
+
+		FREE(dZcoeff);
+	}
+	if(lastIdx) *lastIdx = maxIdx;
+	FREE(icopy);
+	return Zidxs;
+}

zernike_annular.c

@@ -0,0 +1,388 @@
+/*
+ * zernike_annular.c - annular Zernike polynomials
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * 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 2 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, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+
+/*
+ * These polynomials realized according to formulae in:
+@ARTICLE{1981JOSA...71...75M,
+   author = {{Mahajan}, V.~N.},
+    title = "{Zernike annular polynomials for imaging systems with annular pupils.}",
+  journal = {Journal of the Optical Society of America (1917-1983)},
+ keywords = {Astronomical Optics},
+     year = 1981,
+   volume = 71,
+    pages = {75-85},
+   adsurl = {http://adsabs.harvard.edu/abs/1981JOSA...71...75M},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@ARTICLE{2006ApOpt..45.3442H,
+   author = {{Hou}, X. and {Wu}, F. and {Yang}, L. and {Wu}, S. and {Chen}, Q.
+	},
+    title = "{Full-aperture wavefront reconstruction from annular subaperture interferometric data by use of Zernike annular polynomials and a matrix method for testing large aspheric surfaces}",
+  journal = {\ao},
+ keywords = {Mathematics, Optical data processing, Metrology, Optical testing},
+     year = 2006,
+    month = may,
+   volume = 45,
+    pages = {3442-3455},
+      doi = {10.1364/AO.45.003442},
+   adsurl = {http://adsabs.harvard.edu/abs/2006ApOpt..45.3442H},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+ */
+
+#include "zernike.h"
+#include "zern_private.h"
+#include "usefull_macros.h"
+
+// epsilon (==rmin)
+double eps = -1., eps2 = -1.;  // this global variable changes with each conv_r call!
+
+/**
+ * convert coordinates of points on ring into normalized coordinates in [0, 1]
+ * for further calculations of annular Zernike (expression for R_k^0)
+ * @param r0  (i) - original coordinates array
+ * @param Sz      - size of r0
+ * @return dinamically allocated array with new coordinates, which zero element is
+ * 			(0,0)
+ */
+polar *conv_r(polar *r0, int Sz){
+	int x;
+	polar *Pn = MALLOC(polar, Sz);
+	double rmax = -1., rmin = 2.;
+	polar *p = r0, *p1 = Pn;
+	for(x = 0; x < Sz; x++, p++){
+		if(p->r < rmin) rmin = p->r;
+		if(p->r > rmax) rmax = p->r;
+	}
+	if(rmin > rmax || fabs(rmin - rmax) <= DBL_EPSILON || rmax > 1. || rmin < 0.)
+		ERRX(_("polar coordinates should be in [0,1]!"));
+	eps = rmin;
+	eps2 = eps*eps;
+	p = r0;
+	//rmax = 1. - eps2; // 1 - eps^2   -- denominator
+	rmax = rmax*rmax - eps2;
+	for(x = 0; x < Sz; x++, p++, p1++){
+		p1->theta = p->theta;
+		p1->r = sqrt((p->r*p->r - eps2) / rmax); // r = (r^2 - eps^2) / (1 - eps^2)
+	}
+	return Pn;
+}
+
+/**
+ * Allocate 2D array (H arrays of length W)
+ * @param len - array's width  (number of elements in each 1D array
+ * @param num - array's height (number of 1D arrays)
+ * @return dinamically allocated 2D array filled with zeros
+ */
+_2D *_2ALLOC(int len, int num){
+	_2D *r = MALLOC(_2D, 1);
+	r->data = MALLOC(double*, num);
+	int x;
+	for(x = 0; x < num; x++)
+		r->data[x] = MALLOC(double, len);
+	r->len = (size_t)len;
+	r->num = (size_t)num;
+	return r;
+}
+
+void _2FREE(_2D **a){
+	int x, W = (*a)->num;
+	for(x = 0; x < W; x++) free((*a)->data[x]);
+	free((*a)->data);
+	FREE(*a);
+}
+
+/**
+ * Zernike radial function R_n^0
+ * @param n     - order of polynomial
+ * @param Sz    - number of points
+ * @param P (i) - modified polar coordinates: sqrt((r^2-eps^2)/(1-eps^2))
+ * @param R (o) - array with size Sz to which polynomial will be stored
+ */
+void Rn0(int n, int Sz, polar *P, double *R){
+	if(Sz < 1 || !P)
+		ERRX(_("Size of matrix must be > 0!"));
+	if(n < 0)
+		ERRX(_("Order of Zernike polynomial must be > 0!"));
+	if(!FK) build_factorial();
+	int j, k, iup = n/2;
+	double **Rpow = build_rpowR(n, Sz, P);
+	// now fill output matrix
+	double *Zptr = R;
+	polar *p = P;
+	for(j = 0; j < Sz; j++, p++, Zptr++){
+		double Z = 0.;
+		if(p->r > 1.) continue; // throw out points with R>1
+		// calculate R_n^m
+		for(k = 0; k <= iup; k++){ // Sum
+			double F = FK[iup-k];
+			double z = (1. - 2. * (k % 2)) * FK[n - k]    //   (-1)^k * (n-k)!
+				/(//-------------------------------------   ---------------------
+					FK[k]*F*F                             //   k!((n/2-k)!)^2
+				);
+			Z += z * Rpow[n-2*k][j];  // *R^{n-2k}
+		}
+		*Zptr = Z;
+//DBG("order %d, point (%g, %g): %g\n", n, p->r, p->theta, Z);
+	}
+	// free unneeded memory
+	free_rpow(&Rpow, n);
+}
+
+/**
+ * Convert parameters n, m into index in Noll notation, p
+ * @param n, m - Zernike orders
+ * @return order in Noll notation
+ */
+int p_from_nm(int n, int m){
+	return (n*n + 2*n + m)/2;
+}
+
+/**
+ * Create array [p][pt] of annular Zernike polynomials
+ * @param pmax   - max number of polynomial in Noll notation
+ * @param Sz     - size of points array
+ * @param P (i)  - array with points' coordinates
+ * @param Norm(o)- array [0..pmax] of sum^2 for each polynomial (dynamically allocated) or NULL
+ * @return
+ */
+//#undef DBG
+//#define DBG(...)
+_2D *ann_Z(int pmax, int Sz, polar *P, double **Norm){
+	int m, n, j, nmax, mmax, jmax, mW, jm;
+	convert_Zidx(pmax, &nmax, &mmax);
+	mmax = nmax; // for a case when p doesn't include a whole n powers
+	jmax = nmax / 2;
+	mW = mmax + 1; // width of row
+	jm = jmax + mW; // height of Q0 and h
+	polar *Pn = conv_r(P, Sz);
+	if(eps < 0. || eps2 < 0. || eps > 1. || eps2 > 1.)
+		ERRX(_("Wrong epsilon value!!!"));
+	// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+	_2D *h = _2ALLOC(mW, jm); // constants -> array [(jmax+1)+((mmax+1)-1)][mmax+1]
+	_2D *Q0= _2ALLOC(mW, jm); // Q(0) -> array [(jmax+1)+((mmax+1)-1)][mmax+1]
+	_2D *Q = _2ALLOC(Sz, (jmax + 1)*mW); // functions -> array [(jmax+1)*(mmax+1)][Sz]
+DBG("Allocated Q&h, size of Q: %zdx%zd, size of H: %zdx%zd\n", Q->num, Q->len, h->num, h->len);
+DBG("pmax = %d, nmax = mmax = %d, jmax = %d, jm = %d\n", pmax, nmax, jmax, jm);
+DBG("eps=%g, eps2=%g\n", eps, eps2);
+	// in Q index of [j][m][point] is [j*mmax+m][point]
+	double **hd = h->data, **Qd = Q->data, **Q0d = Q0->data;
+	// fill h_j^0 & Q_j^0
+	for(j = 0; j < jm; j++){ // main matrix + triangle
+DBG("j = %d, m = 0\n", j);
+		polar Zero = {0., 0.};
+		// h_j^0 = (1-eps^2)/(2(2j+1))
+		hd[j][0] = (1.-eps2)/(2.*j+1.)/2.;
+		// Q_j^0(rho, eps) = R_{2j}^0(r)
+DBG("Q,h [%d][0], h=%g\n", j, hd[j][0]);
+		if(j <= jmax) Rn0(2*j, Sz, Pn, Qd[mW*j]);
+		Rn0(2*j, 1, &Zero, &Q0d[j][0]); // fill also Q(0)
+	}
+	FREE(Pn);
+	// fill h_j^m & Q_j&m
+	int JMAX = jm-1; // max number of j
+	for(m = 1; m <= mmax; m++, JMAX--){ // m in outern cycle because of Q_{j+1}!
+DBG("\n\nm=%d\n",m);
+		int mminusone = m - 1;
+		int idx = mminusone;// index Q_j^{m-1}
+		for(j = 0; j < JMAX; j++, idx+=mW){
+DBG("\nj=%d\n",j);
+			//             2(2j+2m-1) * h_j^{m-1}
+			// H_j^m = ------------------------------
+			//         (j+m)(1-eps^2) * Q_j^{m-1}(0)
+			double H = 2*(2*(j+m)-1)/(j+m)/(1-eps2) * hd[j][mminusone]/Q0d[j][mminusone];
+			// h_j^m = -H_j^m * Q_{j+1}^{m-1}(0)
+			hd[j][m] = -H * Q0d[j+1][mminusone];
+DBG("hd[%d][%d] = %g (H = %g)\n", j, m, hd[j][m], H);
+			//                     j   Q_i^{m-1}(0) * Q_i^{m-1}(r)
+			// Q_j^m(r) = H_j^m * Sum -----------------------------
+			//                    i=0          h_i^{m-1}
+			if(j <= jmax){
+				int pt;
+				for(pt = 0; pt < Sz; pt++){ // cycle by points
+					int i, idx1 = mminusone;
+					double S = 0;
+					for(i = 0; i <=j; i++, idx1+=mW){ // cycle Sum
+						S += Q0d[i][mminusone] * Qd[idx1][pt] / hd[i][mminusone];
+					}
+					Qd[idx+1][pt] = H * S; // Q_j^m(rho^2)
+//DBG("Q[%d][%d](%d) = %g\n", j, m, pt, H * S);
+				}
+			}
+			// and Q(0):
+			double S = 0.;
+			int i;
+			for(i = 0; i <=j; i++){
+				double qim = Q0d[i][mminusone];
+				S += qim*qim / hd[i][mminusone];
+			}
+			Q0d[j][m] = H * S; // Q_j^m(0)
+			DBG("Q[%d][%d](0) = %g\n", j, m, H * S);
+		}
+	}
+	_2FREE(&Q0);
+
+	// allocate memory for Z_p
+	_2D *Zarr = _2ALLOC(Sz, pmax + 1); // pmax is max number _including_
+	double **Zd = Zarr->data;
+	// now we should calculate R_n^m
+DBG("R_{2j}^0\n");
+	// R_{2j}^0(rho) = Q_j^0(rho^2)
+	size_t S = Sz*sizeof(double);
+	for(j = 0; j <= jmax; j++){
+		int pidx = p_from_nm(2*j, 0);
+		if(pidx > pmax) break;
+		DBG("R[%d]\n", pidx);
+		memcpy(Zd[pidx], Qd[j*mW], S);
+	}
+DBG("R_n^n\n");
+	// R_n^n(rho) = rho^n / sqrt(Sum_{i=0}^n[eps^{2i}])
+	double **Rpow = build_rpowR(nmax, Sz, P); // powers of rho
+	double epssum = 1., epspow = 1.;
+	for(n = 1; n <= nmax; n++){
+		int pidx = p_from_nm(n, -n); // R_n^{-n}
+		if(pidx > pmax) break;
+		DBG("R[%d] (%d, %d)\n", pidx, n, -n);
+		epspow *= eps2; // eps^{2n}
+		epssum += epspow;// sum_0^n eps^{2n}
+		double Sq = sqrt(epssum); // sqrt(sum...)
+		double *rptr = Zd[pidx], *pptr = Rpow[n];
+		int pt;
+		for(pt = 0; pt < Sz; pt++, rptr++, pptr++)
+			*rptr = *pptr / Sq; // rho^n / sqrt(...)
+		int pidx1 = p_from_nm(n, n); // R_n^{n}
+		if(pidx1 > pmax) break;
+		DBG("R[%d] (%d, %d)\n", pidx1, n, n);
+		memcpy(Zd[pidx1], Zd[pidx], S);
+	}
+DBG("R_n^m\n");
+	/*                                   /    1 - eps^2      \
+	   R_{2j+|m|}^{|m|}(rho, eps) = sqrt| ------------------ | * rho^m * Q_j^m(rho^2)
+	                                     \ 2(2j+|m|+1)*h_j^m /                          */
+	for(j = 1; j <= jmax; j++){
+		for(m = 1; m <= mmax; m++){
+			int N = 2*j + m, pt, pidx = p_from_nm(N, -m);
+			if(pidx > pmax) continue;
+			DBG("R[%d]\n", pidx);
+			double *rptr = Zd[pidx], *pptr = Rpow[m], *qptr = Qd[j*mW+m];
+			double hjm = hd[j][m];
+			for(pt = 0; pt < Sz; pt++, rptr++, pptr++, qptr++)
+				*rptr = sqrt((1-eps2)/2/(N+1)/hjm) * (*pptr) * (*qptr);
+			int pidx1 = p_from_nm(N, m);
+			if(pidx1 > pmax) continue;
+			DBG("R[%d]\n", pidx1);
+			memcpy(Zd[pidx1], Zd[pidx], S);
+		}
+	}
+	free_rpow(&Rpow, nmax);
+	_2FREE(&h);
+	_2FREE(&Q);
+	// last step: fill Z_n^m and find sum^2
+	int p;
+	double *sum2 = MALLOC(double, pmax+1);
+	for(p = 0; p <= pmax; p++){
+		convert_Zidx(p, &n, &m);
+		DBG("p=%d (n=%d, m=%d)\n", p, n, m);
+		int pt;
+		double (*fn)(double val);
+		double empty(double _U_ val){return 1.;}
+		if(m == 0) fn = empty;
+		else{if(m > 0)
+			fn = cos;
+		else
+			fn = sin;}
+		double m_abs = fabs((double)m), *rptr = Zd[p], *sptr = &sum2[p];
+		polar *Pptr = P;
+		for(pt = 0; pt < Sz; pt++, rptr++, Pptr++){
+			*rptr *= fn(m_abs * Pptr->theta);
+			//DBG("\tpt %d, val = %g\n", pt, *rptr);
+			*sptr += (*rptr) * (*rptr);
+		}
+		DBG("SUM p^2 = %g\n", *sptr);
+	}
+	if(Norm) *Norm = sum2;
+	else FREE(sum2);
+	return Zarr;
+}
+//#undef DBG
+//#define DBG(...) do{fprintf(stderr, __VA_ARGS__); }while(0)
+
+/**
+ * Compose wavefront by koefficients of annular Zernike polynomials
+ * @params like for ZcomposeR
+ * @return composed wavefront (dynamically allocated)
+ */
+double *ann_Zcompose(int Zsz, double *Zidxs, int Sz, polar *P){
+	int i;
+	double *image = MALLOC(double, Sz);
+	_2D *Zannular = ann_Z(Zsz-1, Sz, P, NULL);
+	double **Zcoeff = Zannular->data;
+	for(i = 0; i < Zsz; i++){ // now we fill array
+		double K = Zidxs[i];
+		if(fabs(K) < Z_prec) continue;
+		int j;
+		double *iptr = image, *zptr = Zcoeff[i];
+		for(j = 0; j < Sz; j++, iptr++, zptr++){
+			*iptr += K * (*zptr); // add next Zernike polynomial
+		}
+	}
+	_2FREE(&Zannular);
+	return image;
+}
+
+
+/**
+ * Decompose wavefront by annular Zernike
+ * @params like ZdecomposeR
+ * @return array of Zernike coefficients
+ */
+double *ann_Zdecompose(int Nmax, int Sz, polar *P, double *heights, int *Zsz, int *lastIdx){
+	int i, pmax, maxIdx = 0;
+	double *Zidxs = NULL, *icopy = NULL;
+	pmax = (Nmax + 1) * (Nmax + 2) / 2; // max Z number in Noll notation
+	Zidxs = MALLOC(double, pmax);
+	icopy = MALLOC(double, Sz);
+	memcpy(icopy, heights, Sz*sizeof(double)); // make image copy to leave it unchanged
+	*Zsz = pmax;
+	double *norm;
+	_2D *Zannular = ann_Z(pmax-1, Sz, P, &norm);
+	double **Zcoeff = Zannular->data;
+	for(i = 0; i < pmax; i++){ // now we fill array
+		int j;
+		double *iptr = icopy, *zptr = Zcoeff[i], K = 0., norm_i = norm[i];
+		for(j = 0; j < Sz; j++, iptr++, zptr++)
+			K += (*zptr) * (*iptr) / norm_i; // multiply matrixes to get coefficient
+		if(fabs(K) < Z_prec)
+			continue; // there's no need to substract values that are less than our precision
+		Zidxs[i] = K;
+		maxIdx = i;
+		iptr = icopy; zptr = Zcoeff[i];
+		for(j = 0; j < Sz; j++, iptr++, zptr++)
+			*iptr -= K * (*zptr); // subtract composed coefficient to reduce high orders values
+	}
+	if(lastIdx) *lastIdx = maxIdx;
+	_2FREE(&Zannular);
+	FREE(norm);
+	FREE(icopy);
+	return Zidxs;
+}