Z coefficients manipulation, save wavefront into png and FITS

wfs_read/Makefile

@@ -1,24 +1,44 @@
+# run `make DEF=...` to add extra defines
 PROGRAM = readwfs
-LDFLAGS = -lm
+LDFLAGS := -fdata-sections -ffunction-sections -Wl,--gc-sections -Wl,--discard-all
+LDFLAGS = -lm -lpng -lcfitsio
 SRCS = $(wildcard *.c)
+DEFINES = $(DEF) -D_XOPEN_SOURCE=1111
+DEFINES += -DEBUG
+CFLAGS = -O2 -fopenmp  -Wall -Werror -Wextra -std=gnu99 -Wno-trampolines
+OBJDIR = mk
+OBJS = $(addprefix $(OBJDIR)/, $(SRCS:%.c=%.o))
+DEPS := $(OBJS:.o=.d)
 CC = gcc
-DEFINES = -D_XOPEN_SOURCE=1111
+
-CXX = gcc
+all : $(OBJDIR) $(PROGRAM)
-CFLAGS = -Wall -Werror -Wextra $(DEFINES) -std=gnu99
+
-OBJS = $(SRCS:.c=.o)
-all : $(PROGRAM)
 $(PROGRAM) : $(OBJS)
-	$(CC) $(CFLAGS) $(OBJS) $(LDFLAGS) -o $(PROGRAM)
+	@echo -e "\t\tLD $(PROGRAM)"
+	$(CC) $(CFLAGS) $(LDFLAGS) $(OBJS) -o $(PROGRAM)
 
-# some addition dependencies
+$(OBJDIR):
-# %.o: %.c
+	mkdir $(OBJDIR)
-#        $(CC) $(LDFLAGS) $(CFLAGS) $< -o $@
-$(SRCS) : %.c : %.h
-	touch $@
 
-%.h: ;
+$(OBJDIR)/%.o: %.c
+	@echo -e "\t\tCC $<"
+	$(CC) -MD -c $(CFLAGS) $(DEFINES) -o $@ $<
 
 clean:
-	/bin/rm -f *.o *~
+	@echo -e "\t\tCLEAN"
-depend:
+	@rm -f $(OBJS)
-	$(CXX) -MM $(CXX.SRCS)
+	@rm -f $(DEPS)
+	@rmdir $(OBJDIR) 2>/dev/null || true
+
+xclean: clean
+	@rm -f $(PROGRAM)
+
+gentags:
+	CFLAGS="$(CFLAGS) $(DEFINES)" geany -g readwfs.c.tags *[hc] 2>/dev/null
+
+.PHONY: gentags clean xclean
+
+ifneq ($(MAKECMDGOALS),clean)
+-include $(DEPS)
+endif
+

wfs_read/README.tags

@@ -0,0 +1,14 @@
+refresh tags:
+
+geany -g readwfs.c.tags *[hc]
+
+after that:
+Tools --> Load Tags File.
+
+read & show matrix:
+M=dlmread('WFS_06012015030446849.matrix'); M([1:3],:)=[];imagesc(M);axis square;
+
+print contour:
+[c h]=contourf(M);
+clabel(c,h)
+axis square; print -dpng wavefront_contour.png

wfs_read/cmdlnopts.c

@@ -40,6 +40,9 @@ glob_pars const Gdefault = {
     ,.wavelength = DEFAULT_WAVELENGTH   // default wavelength
     ,.zzero = 0                         // amount of Z polynomials to be reset
     ,.rotangle = 0.                     // wavefront rotation angle (rotate matrix to -rotangle after computing)
+    ,.tozero = NULL                     // coefficients to be reset (maybe more than one time)
+    ,.addcoef = NULL                    // constant to be added (format x=c, where x is number, c is additive constant)
+    ,.scale = 1.                        // Zernike coefficients' scaling factor
 };
 
 /*
@@ -58,6 +61,9 @@ myoption cmdlnopts[] = {
     {"wavelength",  NEED_ARG, NULL, 'l',    arg_double, APTR(&G.wavelength),_("default wavelength (in meters, microns or nanometers), 101..9999nm")},
     {"zerofirst",   NEED_ARG, NULL, 'z',    arg_int,    APTR(&G.zzero),     _("amount of first Zernike polynomials to be reset to 0")},
     {"rotangle",    NEED_ARG, NULL, 'r',    arg_double, APTR(&G.rotangle),  _("wavefront rotation angle (degrees, CCW)")},
+    {"zero",        MULT_PAR, NULL, '0',    arg_int,    APTR(&G.tozero),    _("reset given Znumber to 0")},
+    {"addconst",    MULT_PAR, NULL, 'a',    arg_string, APTR(&G.addcoef),   _("add constant to given Znumber (e.g. -a4=10 adds 10 to Znum=4)")},
+    {"scale",       NEED_ARG, NULL, 'S',    arg_double, APTR(&G.scale),     _("Zernike coefficients' scaling factor")},
     end_option
 };
 

wfs_read/cmdlnopts.h

@@ -36,7 +36,10 @@ typedef struct{
     char *wfunits;      // units for wavefront measurement in WF map
     double wavelength;  // default wavelength
     int zzero;          // amount of Z polynomials to be reset
+    int **tozero;       // coefficients to be reset (maybe more than one time)
+    char **addcoef;     // constant to be added (format x=c, where x is number, c is additive constant)
     double rotangle;    // wavefront rotation angle (rotate matrix to -rotangle after computing)
+    double scale;       // Zernike coefficients' scaling factor
 } glob_pars;
 
 

wfs_read/main.c

@@ -64,6 +64,9 @@ void proc_DAT(){
             return;
         }
     }
+    if(GP->tozero) z_set_tozero(GP->tozero);
+    if(GP->addcoef) z_set_addcoef(GP->addcoef);
+
     if(strcmp(GP->wfunits, DEFAULT_WF_UNIT)){ // user ask to change default unit
         if(z_set_wfunit(GP->wfunits)){
             WARNX(_("Bad wavefront unit: %s. Should be one of"), GP->wfunits);
@@ -135,6 +138,7 @@ int main(int argc, char** argv){
     initial_setup();
     GP = parse_args(argc, argv);
     if(!GP->inwfs && !GP->indat) ERRX(_("You should give input file name"));
+    z_set_scale(GP->scale);
     if(GP->inwfs){
         proc_WFS();
     }

wfs_read/saveimg.c

@@ -0,0 +1,318 @@
+/*
+ * 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 "saveimg.h"
+#include <png.h>
+#include <fitsio.h>
+#include <math.h>
+#include "zernike.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 + "." + suffix' if outfile != NULL
+ * 		or as 'output + number + "." + suffix' if outfile == NULL
+ * 		number -- any free number from 1 to 9999
+ *
+ * @param outfile (i) - file name
+ * @param suffix  (i) - file suffix
+ * @return created filename (must bee free'd after using) 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(outfile){
+		if(snprintf(buff, 255, "%s%s", outfile, sfx) < 1){
+			DBG("error");
+			return NULL;
+		}
+		return strdup(buff);
+	}else outfile = "output";
+	int num;
+	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, 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, "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)");
+	// OBJECT  / Object name
+	WRITEKEY(TSTRING, "OBJECT", "wavefront", "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 restoration", "Observation program identifier");
+	// AUTHOR / Author of the program
+	WRITEKEY(TSTRING, "AUTHOR", "Edward V. Emelianov <eddy@sao.ru>", "Author of the program");
+	int zerofirst = z_set_Nzero(-1);    // amount of first coefficients set to zero
+	if(zerofirst > 0) 
+		WRITEKEY(TINT, "ZEROFST", &zerofirst, "Amount of first coefficients reset to 0");
+	char *ounit = z_get_wfunit(); // output units
+	if(ounit)
+		WRITEKEY(TSTRING, "OUTPUNIT", ounit, "Output unit");
+	if(0 == strcmp(ounit, "wavelength")){
+		double wavelength = z_get_wavelength(); // wavelength (in meters)
+		WRITEKEY(TDOUBLE, "WAVELEN", &wavelength, "Reference wavelength, m");
+	}
+	double rotangle = z_get_rotangle()*180./M_PI; // wavefront rotation angle (CCW)
+	WRITEKEY(TDOUBLE, "ROTANGLE", &rotangle, "Input wavefront rotation angle (CCW), degr");
+	double step = z_get_step();
+	WRITEKEY(TDOUBLE, "STEP", &step, "Coordinate step on unitary circle");
+	double scale = z_get_scale();
+	WRITEKEY(TDOUBLE, "ZSCALE", &scale, "Scaling factor from original data");
+	int *tozlist = NULL;
+	int tozsz = z_get_tozero(&tozlist);
+	if(tozsz > 0){
+		for(int i = 0; i < tozsz; ++i){
+			char NM[9];
+			int Nul = 0;
+			snprintf(NM, 9, "ZRST%d", tozlist[i]);
+			WRITEKEY(TINT, NM, &Nul, "Coefficient number being reset to 0");
+		}
+	}
+	coeff *addclist = NULL;
+	int addcsz = z_get_addcoef(&addclist);
+	if(addcsz > 0){
+		for(int i = 0; i < addcsz; ++i){
+			char NM[9];
+			snprintf(NM, 9, "ZADD%d", addclist[i].idx);
+			WRITEKEY(TDOUBLE, NM, &addclist[i].addval, "Additional value to given coefficient number");
+		}
+	}
+	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 (i) - filename prefix or NULL to save to "outXXXX.format"
+ * @param sz       - image size (width and height)
+ * @param data (i) - image data
+ * @return number of saved images
+ */
+int writeimg(char *name, size_t sz, double *data){
+	FNAME();
+	char *filename = NULL;
+	int ret = 0;
+	get_stat(data, sz*sz);
+	filename = createfilename(name, "fits");
+	if(filename){
+		ret = writefits(filename, sz, data);
+		FREE(filename);
+	}
+	filename = createfilename(name, "png");
+	if(filename){
+		ret += writepng(filename, sz, data);
+		FREE(filename);
+	}
+	return ret;
+}

wfs_read/saveimg.h

@@ -0,0 +1,28 @@
+/*
+ * 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__
+
+int writeimg(char *name, size_t sz, double *data);
+
+#endif // __SAVEIMG_H__

wfs_read/zernike.c

@@ -24,13 +24,17 @@
 #endif
 #include <math.h>
 #include <strings.h>
+#include <limits.h> // INT_MAX
 #include "zernike.h"
 #include "usefull_macros.h"
+#include "saveimg.h"
 
 #ifndef iabs
 #define iabs(a)  (((a)<(0)) ? (-a) : (a))
 #endif
 
+// scaling factor
+static double zscale = 1.;
 // coordinate step on a grid
 static double coord_step = DEFAULT_CRD_STEP;
 // default wavelength for wavefront (650nm) in meters
@@ -43,15 +47,98 @@ static double *FK = NULL;
 static char *outpunit = DEFAULT_WF_UNIT;
 // amount of first polynomials to reset
 static int Zern_zero = 0;
+// array with indexes to be reset to 0 (copy of G.tozero)
+static int* tozero = NULL;
+static int tozerosz = 0; // size of array
+
+// list of additional coefficients
+static coeff *addcoefflist = NULL;
+static int addcoefsz = 0; // its size
 // matrix rotation angle (in radians)
 static double rotangle = 0.;
+
+/**
+ * setter & getter for zscale
+ */
+void z_set_scale(double s){zscale = s;}
+double z_get_scale(){return zscale;}
+
 /**
  * Set/get value of Zern_zero
  */
 int z_set_Nzero(int val){
-    if(val > -1) Zern_zero  = val;
+    if(val > 0) Zern_zero  = val;
     return Zern_zero;
 }
+/**
+ * Fill list of coefficients to be reset
+ */
+void z_set_tozero(int **val){
+    if(!val) return;
+    size_t cur = tozerosz;
+    while(val[tozerosz]) ++tozerosz;
+    //printf("2zero: %zd\n", tozerosz);
+    tozero = realloc(tozero, tozerosz*sizeof(int));
+    if(!tozero) ERR("realloc()");
+    while(*val){
+        tozero[cur++] = **val++;
+    }
+    //for(cur = 0; cur < tozerosz; ++cur) printf(":: %d\n", tozero[cur]);
+}
+/**
+ * @return size of `tozero` array
+ * @param idxs (o) - `tozero` itself
+ */
+int z_get_tozero(int **idxs){
+    if(!idxs) return -1;
+    *idxs = tozero;
+    return tozerosz;
+}
+
+// parse additional coefficient record (N=val) and fill data
+void parse_coeff(char *rec, coeff *c){
+    if(!rec) ERRX(_("Empty record"));
+    char *s = strdup(rec);
+    if(!s) ERR("strdup()");
+    char *eq = strchr(s, '=');
+    if(!eq || *rec == '=')
+        ERRX(_("Coefficients' records should be like N=val, where N is Znumber, val is its additional value"));
+    *eq++ = 0;
+    char *endptr;
+    long tmp = strtol(s, &endptr, 0);
+    //printf("tmp=%ld, eptr=%s\n", tmp, endptr);
+    if(endptr == rec || *endptr != '\0' || tmp < 0 || tmp > INT_MAX)
+        ERRX(_("Coefficient's index should be a non-negative integer"));
+    c->idx = (int) tmp;
+    if(!str2double(&c->addval, eq)) 
+        ERRX(_("Wrong coefficient"));
+    free(s);
+}
+/**
+ * Fill list of additional coefficients
+ */
+void z_set_addcoef(char **list){
+    if(!list || !*list) return;
+    size_t cur = addcoefsz;
+    while(list[addcoefsz]) ++addcoefsz;
+    //printf("add coeffs: %zd\n", addcoefsz);
+    addcoefflist = realloc(addcoefflist, addcoefsz*sizeof(coeff));
+    if(!addcoefflist) ERR("realloc()");
+    while(*list){
+        parse_coeff(*list++, &addcoefflist[cur++]);
+    }
+    //for(cur = 0; cur < addcoefsz; ++cur) printf("coeff %d += %g\n", addcoefflist[cur].idx, addcoefflist[cur].addval);
+}
+
+/**
+ * @return length of `addcoeflist` array
+ * @param vals (o) - `addcoeflist` itself
+ */
+int z_get_addcoef(coeff **vals){
+    if(!vals) return -1;
+    *vals = addcoefflist;
+    return addcoefsz;
+}
 
 /**
  * Set default coordinate grid step on an unity circle
@@ -90,13 +177,7 @@ double z_get_wavelength(){
     return wavelength;
 }
 
-// for `const char * const *units` thanks to http://stackoverflow.com/a/3875555/1965803
+static wf_units wfunits[] = {
-typedef struct{
-    double wf_coeff;            // multiplier for wavefront units (in .dat files coefficients are in meters)
-    const char * const *units;  // symbol units' names
-} wf_units;
-
-wf_units wfunits[] = {
     { 1.    , (const char * const []){"meter", "m", NULL}},
     {1e3    , (const char * const []){"millimeter", "mm", NULL}},
     {1e6    , (const char * const []){"micrometer", "um", "u", NULL}},
@@ -129,6 +210,10 @@ int z_set_wfunit(char *U){
     return 1;
 }
 
+char *z_get_wfunit(){
+    return outpunit;
+}
+
 double z_get_wfcoeff(){
     return wf_coeff;
 }
@@ -339,10 +424,21 @@ double *Zcompose(int Zsz, double *Zidxs, polcrds *P){
     if(!P || !P->P || !P->Rpow) return NULL;
     int i, Sz = P->Sz;
     for(i = 0; i < Zern_zero; ++i) Zidxs[i] = 0.;
+    for(i = 0; i < tozerosz; i++){
+        int C = tozero[i];
+        if(C > -1 && C < Zsz) Zidxs[C] = 0.;
+        else WARNX(_("Not zero idx %d (should be from 0 to %d)"), C, Zsz-1);
+    }
+    for(i = 0; i < addcoefsz; i++){
+        int C = addcoefflist[i].idx;
+        if(C > -1 && C < Zsz) Zidxs[C] += addcoefflist[i].addval;
+        else WARNX(_("Not change idx %d (should be from 0 to %d)"), C, Zsz-1);
+    }
     double *image = MALLOC(double, Sz);
     for(i = 0; i < Zsz; i++){ // now we fill array
         double K = Zidxs[i];
-        if(fabs(K) < DBL_EPSILON) continue; // 0.0
+        //printf("Z[%d]=%g\n", i, K);
+        if(fabs(K) < DBL_EPSILON) continue; // 0.0m
         int n, m;
         convert_Zidx(i, &n, &m);
         double *Zcoeff = zernfun(n, m, P, NULL);
@@ -389,16 +485,17 @@ int z_save_wavefront(polcrds *P, double *Z, double *std, char *fprefix){
         sum += pt;
         sum2 += pt*pt;
     }
-    sum2 *= wf_coeff, sum *= wf_coeff, max *= wf_coeff, min *= wf_coeff;
+    double coef = wf_coeff * zscale;
-    fprintf(f, "# Wavefront units: %ss, std by all WF: %g, Scope: %g, max: %g, min: %g\n",
+    sum2 *= coef, sum *= coef, max *= coef, min *= coef;
-        outpunit, sum2/Sz + sum*sum/Sz/Sz, max - min, max, min);
+    fprintf(f, "# Wavefront units: %ss, wavelength: %gnm, std by all WF: %g, Scope: %g, max: %g, min: %g\n",
+        outpunit, wavelength*1e9, sum2/Sz + sum*sum/Sz/Sz, max - min, max, min);
     if(Zern_zero) fprintf(f, "# First %d coefficients were cleared\n", Zern_zero);
     fprintf(f, "# X (-1..1)\tY (-1..1)\tZ \tstd_Z\n");
     for(i = 0, z = Z; i < Sz; ++i, ++p, ++z, ++std){
-        double x, y, s, c, r = p->r, zdat = (*z) * wf_coeff;
+        double x, y, s, c, r = p->r, zdat = (*z) * coef;
         sincos(p->theta - rotangle, &s, &c);
         x = r * c, y = r * s;
-        fprintf(f, "%6.3f\t%6.3f\t%9.3g\t%9.3g\n", x, y, zdat, (*std) * wf_coeff);
+        fprintf(f, "%6.3f\t%6.3f\t%9.3g\t%9.3g\n", x, y, zdat, (*std) * coef);
         WF[p->idx] = zdat;
         //DBG("WF[%d] = %g; x=%.1f, y=%.1f", p->idx, zdat,x,y);
     }
@@ -408,15 +505,17 @@ int z_save_wavefront(polcrds *P, double *Z, double *std, char *fprefix){
     printf("try to save to %s\n", filename);
     f = fopen(filename, "w");
     if(!f) goto returning;
-    fprintf(f, "# Wavefront data\n# Units: %ss\n# Step: %g\n", outpunit, coord_step);
+    fprintf(f, "# Wavefront data\n# Units: %ss, wavelength: %gnm\n# Step: %g\n",
+        outpunit, wavelength*1e9, coord_step);
     int x, y;
     // Invert Y axe to have matrix with right Y direction (towards up)
     for(y = WH-1; y > -1; --y){
         double *wptr = &WF[y*WH];
         for(x = 0; x < WH; ++x, ++wptr)
-            fprintf(f, "%6.3f\t", *wptr);
+            fprintf(f, "%6.3g\t", *wptr);
         fprintf(f, "\n");
     }
+    writeimg(fprefix, WH, WF);
     ret = 0;
 returning:
     FREE(filename);
@@ -433,3 +532,6 @@ void z_set_rotangle(double angle){
     printf("angle: %g\n", rotangle*180/M_PI);
 }
 
+double z_get_rotangle(){
+    return rotangle;
+}

wfs_read/zernike.h

@@ -44,6 +44,21 @@ typedef struct{
     int WH;         // Width/Height of matrix
 } polcrds;
 
+// for `const char * const *units` thanks to http://stackoverflow.com/a/3875555/1965803
+typedef struct{
+    double wf_coeff;            // multiplier for wavefront units (in .dat files coefficients are in meters)
+    const char * const *units;  // symbol units' names
+} wf_units;
+
+// additional components to Zernike coefficients
+typedef struct{
+    int idx;        // Znumber
+    double addval;  // additional value
+} coeff;
+
+void z_set_scale(double s);
+double z_get_scale();
+
 int z_set_step(double step);
 double z_get_step();
 
@@ -51,12 +66,21 @@ int z_set_wavelength(double w);
 double z_get_wavelength();
 
 int z_set_wfunit(char *u);
+char *z_get_wfunit();
+
 double z_get_wfcoeff();
 void z_print_wfunits();
 
 void z_set_rotangle(double angle);
+double z_get_rotangle();
 
-int z_set_Nzero(int val);
+int z_set_Nzero(int val); // setter when val > 0 and getter elsewhere
+
+void z_set_tozero(int **val);
+int z_get_tozero(int **idxs);
+
+void z_set_addcoef(char **list);
+int z_get_addcoef(coeff **vals);
 
 void convert_Zidx(int p, int *N, int *M);