pre-alpha

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
+     .inwfs = NULL                      // input WFS file name
-    ,.indat = NULL                     // input DAT file name
+    ,.indat = NULL                      // input DAT file name
-    ,.outname = "wavefront_coords.dat" // output file name
+    ,.outname = NULL                    // output file name prefix (default: basename of input)
-    ,.step = DEFAULT_CRD_STEP          // coordinate step in wavefront map
+    ,.step = DEFAULT_CRD_STEP           // coordinate step in wavefront map
-    ,.wfunits = DEFAULT_WF_UNIT        // units for wavefront measurement in WF map
+    ,.wfunits = DEFAULT_WF_UNIT         // units for wavefront measurement in WF map
-    ,.wavelength = DEFAULT_WAVELENGTH  // default wavelength
+    ,.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
 };
 

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;
 
 

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(GP->zzero) z_set_Nzero(GP->zzero);
-    if(!zerncoeffs){
+    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");
+    polcrds *crds = gen_coords();
-    for(i = 0; i < Zn; ++i) printf("%4d\t%g\n", i, zerncoeffs[i]);
-
-    int L;
-    polar *crds = gen_coords(&L);
     if(!crds){
         WARNX("malloc()");
         return;
     }
-    printf("%d points\n", L);
+    int Sz = crds->Sz;
-    double *surf = Zcompose(Zn, zerncoeffs, L, crds);
+    printf("%d points\n", Sz);
-    if(z_save_wavefront(L, crds, surf, GP->outname))
+    double *surf = MALLOC(double, Sz), *surf2 = MALLOC(double, Sz);
-        WARN(_("Can't open file %s"), GP->outname);
+    for(i = 0; ; ++i){
-    else
+        printf("image %d         \r",i); fflush(stdout);
-        green(_("Saved to %s\n"),  GP->outname);
+        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;
 }
+

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 pointer to mmaped buffer or NULL
- * @return            dynamically allocated array or NULL in case of error
  */
-double *read_dat_file(char *fname, int *sz){
+datfile *open_dat_file(char *fname){
     if(!fname) return NULL;
-    mmapbuf *dbuf = My_mmap(fname);
+    mmapbuf *buf = My_mmap(fname);
-    if(!dbuf) return NULL;
+    if(!buf) return NULL;
-    char *dat = dbuf->data, *eptr = dat + dbuf->len;
+    char *data = buf->data;
-    if(strncasecmp(dat, "time", 4)){
+    if(strncasecmp(data, "time", 4)){
         WARNX(_("Bad header"));
         return NULL;
     }
-    int rd = 0, skipfst = get_hdrval("piston", dat, eptr), L = Z_REALLOC_STEP;
+    char *eptr = data + buf->len;
-    if(skipfst < 0){
+    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;
 }

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__

wfs_read/zernike.c

@@ -19,8 +19,9 @@
  * MA 02110-1301, USA.
  */
 
-
+#ifndef _GNU_SOURCE
 #define _GNU_SOURCE  (1)  // for math.h
+#endif
 #include <math.h>
 #include <strings.h>
 #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;
+    int i, N = n+1;
-	for(i = 0; i < N; i++) FREE((*Rpow)[i]);
+    for(i = 0; i < N; i++) FREE((*Rpow)[i]);
-	FREE(*Rpow);
+    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){
+double *zernfun(int n, int m, polcrds *P, double *norm){
-    if(check_parameters(n, m, Sz, P)) return NULL;
+    if(check_parameters(n, m, P)) return NULL;
-    int j, k, m_abs = iabs(m), iup = (n-m_abs)/2;
+    int j, k, m_abs = iabs(m), iup = (n-m_abs)/2, Sz = P->Sz;
-    double **Rpow = build_rpow(n, Sz, P);
+    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){
+double *Zcompose(int Zsz, double *Zidxs, polcrds *P){
-    int i;
+    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 & R powers
- * @param P (i) - points coordinates
+ * @param Z         (i) - wavefront shift (in lambdas)
- * @param Z (i) - wavefront shift (in lambdas)
+ * @param std       (i) - std of shift in each point
- * @param filename (i) - name of output file
+ * @param filename  (i) - name of output file
  * @return 1 if failed
  */
-int z_save_wavefront(int Sz, polar *P, double *Z, char *filename){
+int z_save_wavefront(polcrds *P, double *Z, double *std, char *fprefix){
-    if(!P || !Z || Sz < 0 || !filename) return 1;
+    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;
+    if(!f) goto returning;
-    fprintf(f, "# X (-1..1)\tY (-1..1)\tZ (%ss)\n", outpunit);
+    int WH = P->WH, max_sz = WH * WH;
-    int i;
+    double *WF = calloc(max_sz, sizeof(double));
-    for(i = 0; i < Sz; ++i, ++P, ++Z){
+    if(!WF) goto returning;
-        double x, y, s, c, r = P->r;
+    // calculate std & scope by all wavefront
-        sincos(P->theta, &s, &c);
+    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);
+}
+

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 *P;       // polar coordinates inside unitary circle
-} polar;
+    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__