diff --git a/wfs_read/cmdlnopts.c b/wfs_read/cmdlnopts.c index 23011af..39539dc 100644 --- a/wfs_read/cmdlnopts.c +++ b/wfs_read/cmdlnopts.c @@ -32,12 +32,14 @@ static glob_pars G; // DEFAULTS // default global parameters glob_pars const Gdefault = { - .inwfs = NULL // input WFS file name - ,.indat = NULL // input DAT file name - ,.outname = "wavefront_coords.dat" // output file name - ,.step = DEFAULT_CRD_STEP // coordinate step in wavefront map - ,.wfunits = DEFAULT_WF_UNIT // units for wavefront measurement in WF map - ,.wavelength = DEFAULT_WAVELENGTH // default wavelength + .inwfs = NULL // input WFS file name + ,.indat = NULL // input DAT file name + ,.outname = NULL // output file name prefix (default: basename of input) + ,.step = DEFAULT_CRD_STEP // coordinate step in wavefront map + ,.wfunits = DEFAULT_WF_UNIT // units for wavefront measurement in WF map + ,.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) }; /* @@ -50,10 +52,12 @@ myoption cmdlnopts[] = { // simple integer parameter with obligatory arg: {"wfs", NEED_ARG, NULL, 'w', arg_string, APTR(&G.inwfs), _("input WFS file name")}, {"dat", NEED_ARG, NULL, 'd', arg_string, APTR(&G.indat), _("input DAT file name")}, - {"output", NEED_ARG, NULL, 'o', arg_string, APTR(&G.outname), _("output file name")}, + {"output", NEED_ARG, NULL, 'o', arg_string, APTR(&G.outname), _("output file name prefix")}, {"step", NEED_ARG, NULL, 's', arg_double, APTR(&G.step), _("coordinate step in wavefront map (R=1)")}, {"wfunits", NEED_ARG, NULL, 'u', arg_string, APTR(&G.wfunits), _("units for wavefront measurement in output WF map")}, {"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)")}, end_option }; diff --git a/wfs_read/cmdlnopts.h b/wfs_read/cmdlnopts.h index 2cf33dc..5d0a76d 100644 --- a/wfs_read/cmdlnopts.h +++ b/wfs_read/cmdlnopts.h @@ -31,10 +31,12 @@ typedef struct{ char *inwfs; // input WFS file name char *indat; // input DAT file name - char *outname; // output file name + char *outname; // output file name prefix double step; // coordinate step in wavefront map char *wfunits; // units for wavefront measurement in WF map double wavelength; // default wavelength + int zzero; // amount of Z polynomials to be reset + double rotangle; // wavefront rotation angle (rotate matrix to -rotangle after computing) } glob_pars; diff --git a/wfs_read/main.c b/wfs_read/main.c index 7d578f7..2917f82 100644 --- a/wfs_read/main.c +++ b/wfs_read/main.c @@ -49,12 +49,14 @@ void proc_WFS(){ } void proc_DAT(){ - int i, Zn; + int i, j, Zn; if(fabs(GP->step - DEFAULT_CRD_STEP) > DBL_EPSILON){ // user change default step - if((i = z_set_step(GP->step))) + if((i = z_set_step(GP->step))){ WARNX(_("Can't change step to %g, value is too %s"), GP->step, i < 0 ? "small" : "big"); return; + } } + if(fabs(GP->rotangle) > DBL_EPSILON) z_set_rotangle(GP->rotangle); if(fabs(GP->wavelength - DEFAULT_WAVELENGTH) > DBL_EPSILON){ // user want to change wavelength // WARNING! This option test should be before changing unit because units depends on wavelength if(z_set_wavelength(GP->wavelength)){ @@ -69,28 +71,61 @@ void proc_DAT(){ return; } } - double *zerncoeffs = read_dat_file(GP->indat, &Zn); - if(!zerncoeffs){ + if(GP->zzero) z_set_Nzero(GP->zzero); + datfile *dat = open_dat_file(GP->indat); + char *fprefix = NULL; + if(!GP->outname){ // default filename + fprefix = strdup(GP->indat); + char *pt = strrchr(fprefix, '.'); + if(pt && pt != fprefix) *pt = 0; + }else fprefix = strdup(GP->outname); + if(!dat){ WARNX(_("Bad DAT file %s"), GP->indat); return; } - printf("Read coefficients:\n"); - for(i = 0; i < Zn; ++i) printf("%4d\t%g\n", i, zerncoeffs[i]); - - int L; - polar *crds = gen_coords(&L); + polcrds *crds = gen_coords(); if(!crds){ WARNX("malloc()"); return; } - printf("%d points\n", L); - double *surf = Zcompose(Zn, zerncoeffs, L, crds); - if(z_save_wavefront(L, crds, surf, GP->outname)) - WARN(_("Can't open file %s"), GP->outname); - else - green(_("Saved to %s\n"), GP->outname); + int Sz = crds->Sz; + printf("%d points\n", Sz); + double *surf = MALLOC(double, Sz), *surf2 = MALLOC(double, Sz); + for(i = 0; ; ++i){ + printf("image %d \r",i); fflush(stdout); + double *zerncoeffs = dat_read_next_line(dat, &Zn); + if(!zerncoeffs) break; // EOF + #ifdef EBUG + //printf("Read coefficients:\n"); + //for(j = 0; j < Zn; ++j) printf("%4d\t%g\n", j, zerncoeffs[j]); + #endif + double *cur = Zcompose(Zn, zerncoeffs, crds); + for(j = 0; j < Sz; ++j){ + double c = cur[j]; + surf[j] += c; + surf2[j] += c*c; + } + free(zerncoeffs); + free(cur); + } + green(_("Got %d iterations, now save file\n"), i); + if(i > 0){ + for(j = 0; j < Sz; ++j){ + surf[j] /= i; // mean + surf2[j] = surf2[j]/i - surf[j]*surf[j]; // std + } + if(z_save_wavefront(crds, surf, surf2, fprefix)) + WARN(_("Can't save files %s"), fprefix); + else + green(_("Saved to %s\n"), fprefix); + }else{ + WARNX(_("Couldn't read any data")); + } + FREE(fprefix); FREE(crds); FREE(surf); + FREE(surf2); + close_dat_file(dat); } /** @@ -108,3 +143,4 @@ int main(int argc, char** argv){ } return 0; } + diff --git a/wfs_read/readdat.c b/wfs_read/readdat.c index 123bf84..69fa618 100644 --- a/wfs_read/readdat.c +++ b/wfs_read/readdat.c @@ -99,32 +99,56 @@ int get_hdrval(char *val, char *dat, char *end){ } /** - * Read .dat file and get Zernike coefficients from it + * Open .dat file * @param fname (i) - .dat file name - * @param sz (o) - size of coefficients' array - * @return dynamically allocated array or NULL in case of error + * @return pointer to mmaped buffer or NULL */ -double *read_dat_file(char *fname, int *sz){ +datfile *open_dat_file(char *fname){ if(!fname) return NULL; - mmapbuf *dbuf = My_mmap(fname); - if(!dbuf) return NULL; - char *dat = dbuf->data, *eptr = dat + dbuf->len; - if(strncasecmp(dat, "time", 4)){ + mmapbuf *buf = My_mmap(fname); + if(!buf) return NULL; + char *data = buf->data; + if(strncasecmp(data, "time", 4)){ WARNX(_("Bad header")); return NULL; } - int rd = 0, skipfst = get_hdrval("piston", dat, eptr), L = Z_REALLOC_STEP; - if(skipfst < 0){ + char *eptr = data + buf->len; + int frst = get_hdrval("piston", data, eptr); + if(frst < 0){ WARNX(_("Dat file don't have OSA Zernike coefficients")); return NULL; } - dat = nextline(dat, eptr); + datfile *dat = MALLOC(datfile, 1); + dat->buf = buf; + dat->eptr = eptr; + dat->curptr = data; + dat->firstcolumn = frst; + return dat; +} + +void close_dat_file(datfile *dat){ + My_munmap(dat->buf); + FREE(dat); +} + +/** + * Read next line from .dat file and get Zernike coefficients from it + * @param dat (i) - input .dat file + * @param sz (o) - size of coefficients' array + * @return dynamically allocated array or NULL in case of error + */ +double *dat_read_next_line(datfile *dat, int *sz){ + if(!dat || !dat->buf){DBG("NULL"); return NULL;} + int rd = 0, L = Z_REALLOC_STEP, skipfst = dat->firstcolumn; + char *buf = dat->curptr, *eptr = dat->eptr; + buf = nextline(buf, eptr); + if(!buf){DBG("EOF"); return NULL;} // EOF double *zern = MALLOC(double, Z_REALLOC_STEP); - while(dat < eptr){ + while(buf < eptr){ double d; - char *next = read_double(dat, eptr, &d); + char *next = read_double(buf, eptr, &d); if(!next) break; - dat = next; + buf = next; if(skipfst > 0){ // skip this value --skipfst; continue; @@ -142,5 +166,6 @@ double *read_dat_file(char *fname, int *sz){ if(*next == '\n') break; } if(sz) *sz = rd; + dat->curptr = buf; return zern; } diff --git a/wfs_read/readdat.h b/wfs_read/readdat.h index 14f300d..e31babc 100644 --- a/wfs_read/readdat.h +++ b/wfs_read/readdat.h @@ -25,6 +25,16 @@ // allocate memory with quantum of this #define Z_REALLOC_STEP (10) -double *read_dat_file(char *fname, int *sz); +typedef struct{ + mmapbuf *buf; // mmaped buffer + char *curptr; // pointer to current symbol + char *eptr; // pointer to end of file + int firstcolumn;// first column with Zernike coefficients + double **Rpow; // powers of R +} datfile; + +double *dat_read_next_line(datfile *dat, int *sz); +datfile *open_dat_file(char *fname); +void close_dat_file(datfile *dat); #endif // __READDAT_H__ diff --git a/wfs_read/zernike.c b/wfs_read/zernike.c index 68016d1..ca96ff1 100644 --- a/wfs_read/zernike.c +++ b/wfs_read/zernike.c @@ -19,8 +19,9 @@ * MA 02110-1301, USA. */ - +#ifndef _GNU_SOURCE #define _GNU_SOURCE (1) // for math.h +#endif #include #include #include "zernike.h" @@ -40,6 +41,17 @@ static double wf_coeff = 1.; static double *FK = NULL; // unit for WF measurement static char *outpunit = DEFAULT_WF_UNIT; +// amount of first polynomials to reset +static int Zern_zero = 0; +// matrix rotation angle (in radians) +static double rotangle = 0.; +/** + * Set/get value of Zern_zero + */ +int z_set_Nzero(int val){ + if(val > -1) Zern_zero = val; + return Zern_zero; +} /** * Set default coordinate grid step on an unity circle @@ -47,7 +59,7 @@ static char *outpunit = DEFAULT_WF_UNIT; * @return 0 if all OK, -1 or 1 if `step` bad */ int z_set_step(double step){ - printf("set to %g\n", step); + printf("Set step to %g\n", step); if(step < DBL_EPSILON) return -1; if(step > 1.) return 1; coord_step = step; @@ -154,69 +166,6 @@ void convert_Zidx(int p, int *N, int *M){ if(N) *N = n; } -/** - * Generate polar coordinates for grid [-1..1] by both coordinates - * with default step - * @param len (o) - size of array - * @return array of coordinates - */ -polar *gen_coords(int *len){ - int WH = 1 + (int)(2. / coord_step), max_sz = WH * WH, L = 0; - polar *coordinates = malloc(max_sz * sizeof(polar)), *cptr = coordinates; - if(!cptr) return NULL; - double x, y; - for(y = -1.; y < 1.; y += coord_step){ - for(x = -1.; x < 1.; x += coord_step){ - double R = sqrt(x*x + y*y); - if(R > 1.) continue; - cptr->r = R; - cptr->theta = atan2(y, x); - ++cptr; - ++L; - } - } - printf("%d points outside circle (ratio = %g, ideal = %g)\n", max_sz - L, ((double)L)/max_sz, M_PI/4.); - if(len) *len = L; - return coordinates; -} - -/** - * Build pre-computed array of factorials from 1 to 100 - */ -void build_factorial(){ - double F = 1.; - int i; - if(FK) return; - FK = MALLOC(double, ZERNIKE_MAX_POWER); - FK[0] = 1.; - for(i = 1; i < ZERNIKE_MAX_POWER; i++) - FK[i] = (F *= (double)i); -} - -/** - * Validation check of zernfun parameters - * return 1 in case of error - */ -int check_parameters(int n, int m, int Sz, polar *P){ - if(Sz < 3 || !P){ - WARNX(_("Size of matrix must be > 2!")); - return 1; - } - if(n > ZERNIKE_MAX_POWER){ - WARNX(_("Order of Zernike polynomial must be <= 100!")); - return 1; - } - int erparm = 0; - if(n < 0) erparm = 1; - if(n < iabs(m)) erparm = 1; // |m| must be <= n - if((n - m) % 2) erparm = 1; // n-m must differ by a prod of 2 - if(erparm) - WARNX(_("Wrong parameters of Zernike polynomial (%d, %d)"), n, m); - else - if(!FK) build_factorial(); - return erparm; -} - /** * Build array with R powers (from 0 to n inclusive) * @param n - power of Zernike polinomial (array size = n+1) @@ -245,28 +194,111 @@ double **build_rpow(int n, int Sz, polar *P){ * @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); + int i, N = n+1; + for(i = 0; i < N; i++) FREE((*Rpow)[i]); + FREE(*Rpow); } +/** + * Free array of coordinates + */ +void free_coords(polcrds *p){ + FREE(p->P); + free_rpow(&p->Rpow, p->N); + free(p); +} +/** + * Generate polar coordinates for grid [-1..1] by both coordinates + * with default step. Correct to rotangle + * @return array of coordinates **without any point outside unitary circle** + */ +polcrds *gen_coords(){ + int WH = 1 + (int)(2. / coord_step), max_sz = WH * WH, L = 0, idx = -1; + polar *coordinates = malloc(max_sz * sizeof(polar)), *cptr = coordinates; + if(!cptr) return NULL; + double x, y, cmax = 1. + coord_step/2.; + for(y = -1.; y < cmax; y += coord_step){ + for(x = -1.; x < cmax; x += coord_step){ + double R = sqrt(x*x + y*y); + ++idx; + if(R > 1.) continue; + cptr->r = R; + cptr->theta = atan2(y, x) + rotangle; + cptr->idx = idx; + ++cptr; + ++L; + } + } + DBG("%d points outside circle (ratio = %g, ideal = %g)", max_sz - L, ((double)L)/max_sz, M_PI/4.); + polcrds *crds = MALLOC(polcrds, 1); + crds->P = coordinates; + crds->Sz = L; + crds->N = 40; // start from 40 + crds->Rpow = build_rpow(crds->N, L, coordinates); + crds->WH = WH; + return crds; +} + +/** + * Build pre-computed array of factorials from 1 to 100 + */ +void build_factorial(){ + double F = 1.; + int i; + if(FK) return; + FK = MALLOC(double, ZERNIKE_MAX_POWER); + FK[0] = 1.; + for(i = 1; i < ZERNIKE_MAX_POWER; i++) + FK[i] = (F *= (double)i); +} + +/** + * Validation check of zernfun parameters + * return 1 in case of error + */ +int check_parameters(int n, int m, polcrds *P){ + if(!P || P->Sz < 3 || !P->P || !P->Rpow){ + WARNX(_("Size of matrix must be > 2!")); + return 1; + } + if(n > ZERNIKE_MAX_POWER){ + WARNX(_("Order of Zernike polynomial must be <= %d!"), ZERNIKE_MAX_POWER); + return 1; + } + int erparm = 0; + if(n < 0) erparm = 1; + if(n < iabs(m)) erparm = 1; // |m| must be <= n + if((n - m) % 2) erparm = 1; // n-m must differ by a prod of 2 + if(erparm) + WARNX(_("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 *zernfun(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_rpow(n, Sz, P); +double *zernfun(int n, int m, polcrds *P, double *norm){ + if(check_parameters(n, m, P)) return NULL; + int j, k, m_abs = iabs(m), iup = (n-m_abs)/2, Sz = P->Sz; + if(n > P->N){ + free_rpow(&P->Rpow, P->N); + P->N = n+10; + P->Rpow = build_rpow(P->N, Sz, P->P); + } double ZSum = 0.; // now fill output matrix double *Zarr = MALLOC(double, Sz); // output matrix double *Zptr = Zarr; - polar *p = P; + polar *p = P->P; + double **Rpow = P->Rpow; for(j = 0; j < Sz; j++, p++, Zptr++){ double Z = 0.; if(p->r > 1.) continue; // throw out points with R>1 @@ -293,8 +325,6 @@ double *zernfun(int n, int m, int Sz, polar *P, double *norm){ ZSum += Z*Z; } if(norm) *norm = ZSum; - // free unneeded memory - free_rpow(&Rpow, n); return Zarr; } @@ -302,18 +332,24 @@ double *zernfun(int n, int m, int Sz, polar *P, double *norm){ * 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) + * @param P(i) - points coordinates & R powers * @return restored image */ -double *Zcompose(int Zsz, double *Zidxs, int Sz, polar *P){ - int i; +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.; 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 int n, m; convert_Zidx(i, &n, &m); - double *Zcoeff = zernfun(n, m, Sz, P, NULL); + double *Zcoeff = zernfun(n, m, P, NULL); + if(!Zcoeff){ + WARNX(_("Can't compute coefficients for n=%d, m=%d!"), n,m); + continue; + } int j; double *iptr = image, *zptr = Zcoeff; for(j = 0; j < Sz; j++, iptr++, zptr++) @@ -325,24 +361,75 @@ double *Zcompose(int Zsz, double *Zidxs, int Sz, polar *P){ /** * Save restored wavefront into file `filename` - * @param Sz - size of `P` - * @param P (i) - points coordinates - * @param Z (i) - wavefront shift (in lambdas) - * @param filename (i) - name of output file + * @param P (i) - points coordinates & R powers + * @param Z (i) - wavefront shift (in lambdas) + * @param std (i) - std of shift in each point + * @param filename (i) - name of output file * @return 1 if failed */ -int z_save_wavefront(int Sz, polar *P, double *Z, char *filename){ - if(!P || !Z || Sz < 0 || !filename) return 1; +int z_save_wavefront(polcrds *P, double *Z, double *std, char *fprefix){ + int Sz = P->Sz, i, ret = 1; + polar *p = P->P; + if(!P || !p || !Z || Sz < 0 || !fprefix) return 1; + /*************** Step 1 - save points coordinates table ***************/ + char *filename = MALLOC(char, strlen(fprefix) + 10); + sprintf(filename, "%s.points", fprefix); + printf("try to save to %s\n", filename); FILE *f = fopen(filename, "w"); - if(!f) return 1; - fprintf(f, "# X (-1..1)\tY (-1..1)\tZ (%ss)\n", outpunit); - int i; - for(i = 0; i < Sz; ++i, ++P, ++Z){ - double x, y, s, c, r = P->r; - sincos(P->theta, &s, &c); + if(!f) goto returning; + int WH = P->WH, max_sz = WH * WH; + double *WF = calloc(max_sz, sizeof(double)); + if(!WF) goto returning; + // calculate std & scope by all wavefront + double *z = Z, sum = 0., sum2 = 0., min = 1e12, max = -1e12; + for(i = 0; i < Sz; ++i, ++z){ + double pt = *z; + if(pt < min) min = pt; + if(pt > max) max = pt; + sum += pt; + sum2 += pt*pt; + } + sum2 *= wf_coeff, sum *= wf_coeff, max *= wf_coeff, min *= wf_coeff; + fprintf(f, "# Wavefront units: %ss, std by all WF: %g, Scope: %g, max: %g, min: %g\n", + outpunit, 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; + sincos(p->theta - rotangle, &s, &c); x = r * c, y = r * s; - fprintf(f, "%g\t%g\t%g\n", x, y, (*Z) * wf_coeff); + fprintf(f, "%6.3f\t%6.3f\t%9.3g\t%9.3g\n", x, y, zdat, (*std) * wf_coeff); + WF[p->idx] = zdat; + //DBG("WF[%d] = %g; x=%.1f, y=%.1f", p->idx, zdat,x,y); } fclose(f); - return 0; + /*************** Step 2 - save matrix of data ***************/ + sprintf(filename, "%s.matrix", 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); + 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, "\n"); + } + ret = 0; +returning: + FREE(filename); + return ret; } + + +/** + * change rotation angle (0..2pi) + */ +void z_set_rotangle(double angle){ + angle -= floor(angle/360.) * 360.; + rotangle = angle * M_PI/180.; + printf("angle: %g\n", rotangle*180/M_PI); +} + diff --git a/wfs_read/zernike.h b/wfs_read/zernike.h index 5c924b8..ab9fd94 100644 --- a/wfs_read/zernike.h +++ b/wfs_read/zernike.h @@ -31,10 +31,18 @@ // max power of Zernike polynomial #define ZERNIKE_MAX_POWER (100) +typedef struct{ + double r,theta; // polar coordinates + int idx; // index of given point in square matrix +} polar; typedef struct{ - double r,theta; -} polar; + polar *P; // polar coordinates inside unitary circle + double **Rpow; // powers of R + int N; // max power of Zernike coeffs + int Sz; // size of P + int WH; // Width/Height of matrix +} polcrds; int z_set_step(double step); double z_get_step(); @@ -46,11 +54,18 @@ int z_set_wfunit(char *u); double z_get_wfcoeff(); void z_print_wfunits(); +void z_set_rotangle(double angle); + +int z_set_Nzero(int val); + void convert_Zidx(int p, int *N, int *M); -polar *gen_coords(int *len); +polcrds *gen_coords(); +void free_coords(polcrds *p); -double *Zcompose(int Zsz, double *Zidxs, int Sz, polar *P); +double **build_rpow(int n, int Sz, polar *P); -int z_save_wavefront(int Sz, polar *P, double *Z, char *filename); +double *Zcompose(int Zsz, double *Zidxs, polcrds *P); + +int z_save_wavefront(polcrds *P, double *Z, double *std, char *fprefix); #endif // __ZERNIKE_H__