#ifndef muscatt_functions_h
#define muscatt_functions_h

//[included Functions]
#define fft2d_alloc_error_check(p) { \
    if ((p) == NULL) { \
        fprintf(stderr, "fft2d memory allocation error\n"); \
        exit(1); \
    } \
}

#ifdef USE_FFT2D_PTHREADS
#define USE_FFT2D_THREADS
#ifndef FFT2D_MAX_THREADS
#define FFT2D_MAX_THREADS 4
#endif
#ifndef FFT2D_THREADS_BEGIN_N
#define FFT2D_THREADS_BEGIN_N 65536
#endif
#include <pthread.h>
#define fft2d_thread_t pthread_t
#define fft2d_thread_create(thp,func,argp) { \
    if (pthread_create(thp, NULL, func, (void *) (argp)) != 0) { \
        fprintf(stderr, "fft2d thread error\n"); \
        exit(1); \
    } \
}
#define fft2d_thread_wait(th) { \
    if (pthread_join(th, NULL) != 0) { \
        fprintf(stderr, "fft2d thread error\n"); \
        exit(1); \
    } \
}
#endif /* USE_FFT2D_PTHREADS */

#ifdef USE_FFT2D_WINTHREADS
#define USE_FFT2D_THREADS
#ifndef FFT2D_MAX_THREADS
#define FFT2D_MAX_THREADS 4
#endif
#ifndef FFT2D_THREADS_BEGIN_N
#define FFT2D_THREADS_BEGIN_N 131072
#endif
#include <windows.h>
#define fft2d_thread_t HANDLE
#define fft2d_thread_create(thp,func,argp) { \
    DWORD thid; \
    *(thp) = CreateThread(NULL, 0, (LPTHREAD_START_ROUTINE) (func), (LPVOID) (argp), 0, &thid); \
    if (*(thp) == 0) { \
        fprintf(stderr, "fft2d thread error\n"); \
        exit(1); \
    } \
}
#define fft2d_thread_wait(th) { \
    WaitForSingleObject(th, INFINITE); \
    CloseHandle(th); \
}
#endif /* USE_FFT2D_WINTHREADS */





double min(double A, double B) {
if (A<B) return A; else return B;
};

double max(double A, double B) {
if (A>B) return A; else return B;
};

int imin(int A, int B) {
if (A<B) return A; else return B;
};

int imax(int A, int B) {
if (A>B) return A; else return B;
};

void linreg(double X[100], double Y[100], double Ye[100], int N, double par[3]) {
  double wt,sX,sXX,sY,sYY,sXY,sW;
  int i;
  sX = 0.0;
  sXX= 0.0;
  sY = 0.0;
  sYY= 0.0;
  sXY= 0.0;
  sW = 0.0;
  for (i=0;i<imin(N,12);i++) {
    wt  = 1/(Ye[i]*Ye[i]);
    sX += X[i]*wt;
    sXX+= X[i]*X[i]*wt;
    sY += Y[i]*wt;
    sYY+= Y[i]*Y[i]*wt;
    sXY+= X[i]*Y[i]*wt;
    sW += wt;
  };
  sX  /= sW;
  sXX /= sW;
  sY  /= sW;
  sYY /= sW;
  sXY /= sW;
  par[0] = (sXY-sX*sY)/sqrt((sXX-sX*sX)*(sYY-sY*sY));         /* correlation coefficient */
  par[1] = par[0]*sqrt((sYY-sY*sY)/(sXX-sX*sX));              /* slope */
  par[2] = sY - sX*par[1];                                    /* ordinate value */
  par[0] = fabs(par[0]);
  /* fprintf(stdout,"%i %12.4e %12.4e %12.4e \n",N,par[0],par[1],par[2]); */
};

void radavg(int steps, int nd, int pos, double Qmin, double **b, double *Qc, double In[300][16]) {
  int i,j;
  double pi,rad,phi,xs,ys;
  int phisteps,ix1,iy1,ix2,iy2;
  double wt, twt, avg;
  int anz;
  double rex,rey;

  pi       = 3.141592653589793238462643;      /* constant */

  for (i=0;i<nd;i++) {

    rad = (steps-1)*(Qc[i]-Qmin)/(Qc[nd-1]-Qmin);
    phisteps = imax((int)(-floor(-0.5*pi*rad)+0.5),1);

    twt = 0.0;      /* totales Gewicht    >>>RESET<<<               */
    avg = 0.0;      /* Mittel                                       */
    anz = 0;        /* Anzahl der Punkte (vielleicht nicht n"otig)  */

    for (j=0;j<phisteps;j++) {

      phi = 0.5*pi*(j+0.5)/phisteps;
      xs  = rad * cos(phi);
      ys  = rad * sin(phi);

      ix1 = floor(xs);         /*  Eckpunkte */
      iy1 = floor(ys);
      ix2 = ix1 + 1;
      iy2 = iy1 + 1;
      rex = xs - ix1;          /*  Reste f"ur Gewichtung sp"ater */
      rey = ys - iy1;

      if (ix1>=0 && ix1<steps && iy1>=0 && iy1<steps) {
        wt   = (1.0-rex)*(1.0-rey);
        twt += wt;
        avg += wt*b[ix1][iy1];
        anz++;
      };
      if (ix2>=0 && ix2<steps && iy1>=0 && iy1<steps) {
        wt   = (    rex)*(1.0-rey);
        twt += wt;
        avg += wt*b[ix2][iy1];
        anz++;
      };
      if (ix1>=0 && ix1<steps && iy2>=0 && iy2<steps) {
        wt   = (1.0-rex)*(    rey);
        twt += wt;
        avg += wt*b[ix1][iy2];
        anz++;
      };
      if (ix2>=0 && ix2<steps && iy2>=0 && iy2<steps) {
        wt   = (    rex)*(    rey);
        twt += wt;
        avg += wt*b[ix2][iy2];
        anz++;
      };
    };
    if (twt>0.0) {avg /= twt;};
    In[i][pos]=avg;
  }; 
};

void radseg(int steps, int nd, int k, int pos, double Qmin, double **b, double *Qc, double In[300][16]) {
  int i,j;
  double pi,rad,phi,xs,ys;
  int phisteps,ix1,iy1,ix2,iy2;
  double wt, twt, avg;
  int anz;
  double rex,rey;

  pi       = 3.141592653589793238462643;      /* constant */

  i=k;

    rad = (steps-1)*(Qc[i]-Qmin)/(Qc[nd-1]-Qmin);
    phisteps = imax((int)(-floor(-0.5*pi*rad)+0.5),1);

    twt = 0.0;      /* totales Gewicht    >>>RESET<<<               */
    avg = 0.0;      /* Mittel                                       */
    anz = 0;        /* Anzahl der Punkte (vielleicht nicht n"otig)  */

    for (j=0;j<phisteps;j++) {

      phi = 0.5*pi*(j+0.5)/phisteps;
      xs  = rad * cos(phi);
      ys  = rad * sin(phi);

      ix1 = floor(xs);         /*  Eckpunkte */
      iy1 = floor(ys);
      ix2 = ix1 + 1;
      iy2 = iy1 + 1;
      rex = xs - ix1;          /*  Reste f"ur Gewichtung sp"ater */
      rey = ys - iy1;

      if (ix1>=0 && ix1<steps && iy1>=0 && iy1<steps) {
        wt   = (1.0-rex)*(1.0-rey);
        twt += wt;
        avg += wt*b[ix1][iy1];
        anz++;
      };
      if (ix2>=0 && ix2<steps && iy1>=0 && iy1<steps) {
        wt   = (    rex)*(1.0-rey);
        twt += wt;
        avg += wt*b[ix2][iy1];
        anz++;
      };
      if (ix1>=0 && ix1<steps && iy2>=0 && iy2<steps) {
        wt   = (1.0-rex)*(    rey);
        twt += wt;
        avg += wt*b[ix1][iy2];
        anz++;
      };
      if (ix2>=0 && ix2<steps && iy2>=0 && iy2<steps) {
        wt   = (    rex)*(    rey);
        twt += wt;
        avg += wt*b[ix2][iy2];
        anz++;
      };
    };
    if (twt>0.0) {avg /= twt;};
    In[i][pos]=avg;
};

/* ...................................................................................................................... */

/*
Fast Fourier/Cosine/Sine Transform
    dimension   :two
    data length :power of 2
    decimation  :frequency
    radix       :split-radix, row-column
    data        :inplace
    table       :use
functions
    cdft2d: Complex Discrete Fourier Transform
    rdft2d: Real Discrete Fourier Transform
    ddct2d: Discrete Cosine Transform
    ddst2d: Discrete Sine Transform
function prototypes
    void cdft2d(int, int, int, double **, double *, int *, double *);
    void rdft2d(int, int, int, double **, double *, int *, double *);
    void rdft2dsort(int, int, int, double **);
    void ddct2d(int, int, int, double **, double *, int *, double *);
    void ddst2d(int, int, int, double **, double *, int *, double *);
necessary package
    fftsg.c  : 1D-FFT package
macro definitions
    USE_FFT2D_PTHREADS : default=not defined
        FFT2D_MAX_THREADS     : must be 2^N, default=4
        FFT2D_THREADS_BEGIN_N : default=65536
    USE_FFT2D_WINTHREADS : default=not defined
        FFT2D_MAX_THREADS     : must be 2^N, default=4
        FFT2D_THREADS_BEGIN_N : default=131072


-------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
    [definition]
        <case1> IDCT (excluding scale)
            C[k1][k2] = sum_j1=0^n1-1 sum_j2=0^n2-1 a[j1][j2] * 
                            cos(pi*j1*(k1+1/2)/n1) * 
                            cos(pi*j2*(k2+1/2)/n2), 
                            0<=k1<n1, 0<=k2<n2
        <case2> DCT
            C[k1][k2] = sum_j1=0^n1-1 sum_j2=0^n2-1 a[j1][j2] * 
                            cos(pi*(j1+1/2)*k1/n1) * 
                            cos(pi*(j2+1/2)*k2/n2), 
                            0<=k1<n1, 0<=k2<n2
    [usage]
        <case1>
            ip[0] = 0; // first time only
            ddct2d(n1, n2, 1, a, t, ip, w);
        <case2>
            ip[0] = 0; // first time only
            ddct2d(n1, n2, -1, a, t, ip, w);
    [parameters]
        n1     :data length (int)
                n1 >= 2, n1 = power of 2
        n2     :data length (int)
                n2 >= 2, n2 = power of 2
        a[0...n1-1][0...n2-1]
               :input/output data (double **)
                output data
                    a[k1][k2] = C[k1][k2], 0<=k1<n1, 0<=k2<n2
        t[0...*]
               :work area (double *)
                length of t >= 4*n1,                   if single thread, 
                length of t >= 4*n1*FFT2D_MAX_THREADS, if multi threads, 
                t is dynamically allocated, if t == NULL.
        ip[0...*]
               :work area for bit reversal (int *)
                length of ip >= 2+sqrt(n)
                (n = max(n1/2, n2/2))
                ip[0],ip[1] are pointers of the cos/sin table.
        w[0...*]
               :cos/sin table (double *)
                length of w >= max(n1*3/2, n2*3/2)
                w[],ip[] are initialized if ip[0] == 0.
    [remark]
        Inverse of 
            ddct2d(n1, n2, -1, a, t, ip, w);
        is 
            for (j1 = 0; j1 <= n1 - 1; j1++) {
                a[j1][0] *= 0.5;
            }
            for (j2 = 0; j2 <= n2 - 1; j2++) {
                a[0][j2] *= 0.5;
            }
            ddct2d(n1, n2, 1, a, t, ip, w);
            for (j1 = 0; j1 <= n1 - 1; j1++) {
                for (j2 = 0; j2 <= n2 - 1; j2++) {
                    a[j1][j2] *= 4.0 / n1 / n2;
                }
            }
        .
*/

void ddct2d(int n1, int n2, int isgn, double **a, double *t, 
    int *ip, double *w)
{
    void makewt(int nw, int *ip, double *w);
    void makect(int nc, int *ip, double *c);
    void ddct(int n, int isgn, double *a, int *ip, double *w);
    void ddxt2d_sub(int n1, int n2, int ics, int isgn, double **a, 
        double *t, int *ip, double *w);
#ifdef USE_FFT2D_THREADS
    void ddxt2d0_subth(int n1, int n2, int ics, int isgn, double **a, 
        int *ip, double *w);
    void ddxt2d_subth(int n1, int n2, int ics, int isgn, double **a, 
        double *t, int *ip, double *w);
#endif /* USE_FFT2D_THREADS */
    int n, nw, nc, itnull, nthread, nt, i;
    
    n = n1;
    if (n < n2) {
        n = n2;
    }
    nw = ip[0];
    if (n > (nw << 2)) {
        nw = n >> 2;
        makewt(nw, ip, w);
    }
    nc = ip[1];
    if (n > nc) {
        nc = n;
        makect(nc, ip, w + nw);
    }
    itnull = 0;
    if (t == NULL) {
        itnull = 1;
        nthread = 1;
#ifdef USE_FFT2D_THREADS
        nthread = FFT2D_MAX_THREADS;
#endif /* USE_FFT2D_THREADS */
        nt = 4 * nthread * n1;
        if (n2 == 2 * nthread) {
            nt >>= 1;
        } else if (n2 < 2 * nthread) {
            nt >>= 2;
        }
        t = (double *) malloc(sizeof(double) * nt);
        fft2d_alloc_error_check(t);
    }
#ifdef USE_FFT2D_THREADS
    if ((double) n1 * n2 >= (double) FFT2D_THREADS_BEGIN_N) {
        ddxt2d0_subth(n1, n2, 0, isgn, a, ip, w);
        ddxt2d_subth(n1, n2, 0, isgn, a, t, ip, w);
    } else 
#endif /* USE_FFT2D_THREADS */
    {
        for (i = 0; i < n1; i++) {
            ddct(n2, isgn, a[i], ip, w);
        }
        ddxt2d_sub(n1, n2, 0, isgn, a, t, ip, w);
    }
    if (itnull != 0) {
        free(t);
    }
}



/* -------- child routines -------- */


void cdft2d_sub(int n1, int n2, int isgn, double **a, double *t, 
    int *ip, double *w)
{
    void cdft(int n, int isgn, double *a, int *ip, double *w);
    int i, j;
    
    if (n2 > 4) {
        for (j = 0; j < n2; j += 8) {
            for (i = 0; i < n1; i++) {
                t[2 * i] = a[i][j];
                t[2 * i + 1] = a[i][j + 1];
                t[2 * n1 + 2 * i] = a[i][j + 2];
                t[2 * n1 + 2 * i + 1] = a[i][j + 3];
                t[4 * n1 + 2 * i] = a[i][j + 4];
                t[4 * n1 + 2 * i + 1] = a[i][j + 5];
                t[6 * n1 + 2 * i] = a[i][j + 6];
                t[6 * n1 + 2 * i + 1] = a[i][j + 7];
            }
            cdft(2 * n1, isgn, t, ip, w);
            cdft(2 * n1, isgn, &t[2 * n1], ip, w);
            cdft(2 * n1, isgn, &t[4 * n1], ip, w);
            cdft(2 * n1, isgn, &t[6 * n1], ip, w);
            for (i = 0; i < n1; i++) {
                a[i][j] = t[2 * i];
                a[i][j + 1] = t[2 * i + 1];
                a[i][j + 2] = t[2 * n1 + 2 * i];
                a[i][j + 3] = t[2 * n1 + 2 * i + 1];
                a[i][j + 4] = t[4 * n1 + 2 * i];
                a[i][j + 5] = t[4 * n1 + 2 * i + 1];
                a[i][j + 6] = t[6 * n1 + 2 * i];
                a[i][j + 7] = t[6 * n1 + 2 * i + 1];
            }
        }
    } else if (n2 == 4) {
        for (i = 0; i < n1; i++) {
            t[2 * i] = a[i][0];
            t[2 * i + 1] = a[i][1];
            t[2 * n1 + 2 * i] = a[i][2];
            t[2 * n1 + 2 * i + 1] = a[i][3];
        }
        cdft(2 * n1, isgn, t, ip, w);
        cdft(2 * n1, isgn, &t[2 * n1], ip, w);
        for (i = 0; i < n1; i++) {
            a[i][0] = t[2 * i];
            a[i][1] = t[2 * i + 1];
            a[i][2] = t[2 * n1 + 2 * i];
            a[i][3] = t[2 * n1 + 2 * i + 1];
        }
    } else if (n2 == 2) {
        for (i = 0; i < n1; i++) {
            t[2 * i] = a[i][0];
            t[2 * i + 1] = a[i][1];
        }
        cdft(2 * n1, isgn, t, ip, w);
        for (i = 0; i < n1; i++) {
            a[i][0] = t[2 * i];
            a[i][1] = t[2 * i + 1];
        }
    }
}


void rdft2d_sub(int n1, int n2, int isgn, double **a)
{
    int n1h, i, j;
    double xi;
    
    n1h = n1 >> 1;
    if (isgn < 0) {
        for (i = 1; i < n1h; i++) {
            j = n1 - i;
            xi = a[i][0] - a[j][0];
            a[i][0] += a[j][0];
            a[j][0] = xi;
            xi = a[j][1] - a[i][1];
            a[i][1] += a[j][1];
            a[j][1] = xi;
        }
    } else {
        for (i = 1; i < n1h; i++) {
            j = n1 - i;
            a[j][0] = 0.5 * (a[i][0] - a[j][0]);
            a[i][0] -= a[j][0];
            a[j][1] = 0.5 * (a[i][1] + a[j][1]);
            a[i][1] -= a[j][1];
        }
    }
}


void ddxt2d_sub(int n1, int n2, int ics, int isgn, double **a, 
    double *t, int *ip, double *w)
{
    void ddct(int n, int isgn, double *a, int *ip, double *w);
    void ddst(int n, int isgn, double *a, int *ip, double *w);
    int i, j;
    
    if (n2 > 2) {
        for (j = 0; j < n2; j += 4) {
            for (i = 0; i < n1; i++) {
                t[i] = a[i][j];
                t[n1 + i] = a[i][j + 1];
                t[2 * n1 + i] = a[i][j + 2];
                t[3 * n1 + i] = a[i][j + 3];
            }
            if (ics == 0) {
                ddct(n1, isgn, t, ip, w);
                ddct(n1, isgn, &t[n1], ip, w);
                ddct(n1, isgn, &t[2 * n1], ip, w);
                ddct(n1, isgn, &t[3 * n1], ip, w);
            } else {
                ddst(n1, isgn, t, ip, w);
                ddst(n1, isgn, &t[n1], ip, w);
                ddst(n1, isgn, &t[2 * n1], ip, w);
                ddst(n1, isgn, &t[3 * n1], ip, w);
            }
            for (i = 0; i < n1; i++) {
                a[i][j] = t[i];
                a[i][j + 1] = t[n1 + i];
                a[i][j + 2] = t[2 * n1 + i];
                a[i][j + 3] = t[3 * n1 + i];
            }
        }
    } else if (n2 == 2) {
        for (i = 0; i < n1; i++) {
            t[i] = a[i][0];
            t[n1 + i] = a[i][1];
        }
        if (ics == 0) {
            ddct(n1, isgn, t, ip, w);
            ddct(n1, isgn, &t[n1], ip, w);
        } else {
            ddst(n1, isgn, t, ip, w);
            ddst(n1, isgn, &t[n1], ip, w);
        }
        for (i = 0; i < n1; i++) {
            a[i][0] = t[i];
            a[i][1] = t[n1 + i];
        }
    }
}



#ifdef USE_FFT2D_THREADS
struct fft2d_arg_st {
    int nthread;
    int n0;
    int n1;
    int n2;
    int ic;
    int isgn;
    double **a;
    double *t;
    int *ip;
    double *w;
};
typedef struct fft2d_arg_st fft2d_arg_t;


void xdft2d0_subth(int n1, int n2, int icr, int isgn, double **a, 
    int *ip, double *w)
{
    void *xdft2d0_th(void *p);
    fft2d_thread_t th[FFT2D_MAX_THREADS];
    fft2d_arg_t ag[FFT2D_MAX_THREADS];
    int nthread, i;
    
    nthread = FFT2D_MAX_THREADS;
    if (nthread > n1) {
        nthread = n1;
    }
    for (i = 0; i < nthread; i++) {
        ag[i].nthread = nthread;
        ag[i].n0 = i;
        ag[i].n1 = n1;
        ag[i].n2 = n2;
        ag[i].ic = icr;
        ag[i].isgn = isgn;
        ag[i].a = a;
        ag[i].ip = ip;
        ag[i].w = w;
        fft2d_thread_create(&th[i], xdft2d0_th, &ag[i]);
    }
    for (i = 0; i < nthread; i++) {
        fft2d_thread_wait(th[i]);
    }
}


void cdft2d_subth(int n1, int n2, int isgn, double **a, double *t, 
    int *ip, double *w)
{
    void *cdft2d_th(void *p);
    fft2d_thread_t th[FFT2D_MAX_THREADS];
    fft2d_arg_t ag[FFT2D_MAX_THREADS];
    int nthread, nt, i;
    
    nthread = FFT2D_MAX_THREADS;
    nt = 8 * n1;
    if (n2 == 4 * FFT2D_MAX_THREADS) {
        nt >>= 1;
    } else if (n2 < 4 * FFT2D_MAX_THREADS) {
        nthread = n2 >> 1;
        nt >>= 2;
    }
    for (i = 0; i < nthread; i++) {
        ag[i].nthread = nthread;
        ag[i].n0 = i;
        ag[i].n1 = n1;
        ag[i].n2 = n2;
        ag[i].isgn = isgn;
        ag[i].a = a;
        ag[i].t = &t[nt * i];
        ag[i].ip = ip;
        ag[i].w = w;
        fft2d_thread_create(&th[i], cdft2d_th, &ag[i]);
    }
    for (i = 0; i < nthread; i++) {
        fft2d_thread_wait(th[i]);
    }
}


void ddxt2d0_subth(int n1, int n2, int ics, int isgn, double **a, 
    int *ip, double *w)
{
    void *ddxt2d0_th(void *p);
    fft2d_thread_t th[FFT2D_MAX_THREADS];
    fft2d_arg_t ag[FFT2D_MAX_THREADS];
    int nthread, i;
    
    nthread = FFT2D_MAX_THREADS;
    if (nthread > n1) {
        nthread = n1;
    }
    for (i = 0; i < nthread; i++) {
        ag[i].nthread = nthread;
        ag[i].n0 = i;
        ag[i].n1 = n1;
        ag[i].n2 = n2;
        ag[i].ic = ics;
        ag[i].isgn = isgn;
        ag[i].a = a;
        ag[i].ip = ip;
        ag[i].w = w;
        fft2d_thread_create(&th[i], ddxt2d0_th, &ag[i]);
    }
    for (i = 0; i < nthread; i++) {
        fft2d_thread_wait(th[i]);
    }
}


void ddxt2d_subth(int n1, int n2, int ics, int isgn, double **a, 
    double *t, int *ip, double *w)
{
    void *ddxt2d_th(void *p);
    fft2d_thread_t th[FFT2D_MAX_THREADS];
    fft2d_arg_t ag[FFT2D_MAX_THREADS];
    int nthread, nt, i;
    
    nthread = FFT2D_MAX_THREADS;
    nt = 4 * n1;
    if (n2 == 2 * FFT2D_MAX_THREADS) {
        nt >>= 1;
    } else if (n2 < 2 * FFT2D_MAX_THREADS) {
        nthread = n2;
        nt >>= 2;
    }
    for (i = 0; i < nthread; i++) {
        ag[i].nthread = nthread;
        ag[i].n0 = i;
        ag[i].n1 = n1;
        ag[i].n2 = n2;
        ag[i].ic = ics;
        ag[i].isgn = isgn;
        ag[i].a = a;
        ag[i].t = &t[nt * i];
        ag[i].ip = ip;
        ag[i].w = w;
        fft2d_thread_create(&th[i], ddxt2d_th, &ag[i]);
    }
    for (i = 0; i < nthread; i++) {
        fft2d_thread_wait(th[i]);
    }
}


void *xdft2d0_th(void *p)
{
    void cdft(int n, int isgn, double *a, int *ip, double *w);
    void rdft(int n, int isgn, double *a, int *ip, double *w);
    int nthread, n0, n1, n2, icr, isgn, *ip, i;
    double **a, *w;
    
    nthread = ((fft2d_arg_t *) p)->nthread;
    n0 = ((fft2d_arg_t *) p)->n0;
    n1 = ((fft2d_arg_t *) p)->n1;
    n2 = ((fft2d_arg_t *) p)->n2;
    icr = ((fft2d_arg_t *) p)->ic;
    isgn = ((fft2d_arg_t *) p)->isgn;
    a = ((fft2d_arg_t *) p)->a;
    ip = ((fft2d_arg_t *) p)->ip;
    w = ((fft2d_arg_t *) p)->w;
    if (icr == 0) {
        for (i = n0; i < n1; i += nthread) {
            cdft(n2, isgn, a[i], ip, w);
        }
    } else {
        for (i = n0; i < n1; i += nthread) {
            rdft(n2, isgn, a[i], ip, w);
        }
    }
    return (void *) 0;
}


void *cdft2d_th(void *p)
{
    void cdft(int n, int isgn, double *a, int *ip, double *w);
    int nthread, n0, n1, n2, isgn, *ip, i, j;
    double **a, *t, *w;
    
    nthread = ((fft2d_arg_t *) p)->nthread;
    n0 = ((fft2d_arg_t *) p)->n0;
    n1 = ((fft2d_arg_t *) p)->n1;
    n2 = ((fft2d_arg_t *) p)->n2;
    isgn = ((fft2d_arg_t *) p)->isgn;
    a = ((fft2d_arg_t *) p)->a;
    t = ((fft2d_arg_t *) p)->t;
    ip = ((fft2d_arg_t *) p)->ip;
    w = ((fft2d_arg_t *) p)->w;
    if (n2 > 4 * nthread) {
        for (j = 8 * n0; j < n2; j += 8 * nthread) {
            for (i = 0; i < n1; i++) {
                t[2 * i] = a[i][j];
                t[2 * i + 1] = a[i][j + 1];
                t[2 * n1 + 2 * i] = a[i][j + 2];
                t[2 * n1 + 2 * i + 1] = a[i][j + 3];
                t[4 * n1 + 2 * i] = a[i][j + 4];
                t[4 * n1 + 2 * i + 1] = a[i][j + 5];
                t[6 * n1 + 2 * i] = a[i][j + 6];
                t[6 * n1 + 2 * i + 1] = a[i][j + 7];
            }
            cdft(2 * n1, isgn, t, ip, w);
            cdft(2 * n1, isgn, &t[2 * n1], ip, w);
            cdft(2 * n1, isgn, &t[4 * n1], ip, w);
            cdft(2 * n1, isgn, &t[6 * n1], ip, w);
            for (i = 0; i < n1; i++) {
                a[i][j] = t[2 * i];
                a[i][j + 1] = t[2 * i + 1];
                a[i][j + 2] = t[2 * n1 + 2 * i];
                a[i][j + 3] = t[2 * n1 + 2 * i + 1];
                a[i][j + 4] = t[4 * n1 + 2 * i];
                a[i][j + 5] = t[4 * n1 + 2 * i + 1];
                a[i][j + 6] = t[6 * n1 + 2 * i];
                a[i][j + 7] = t[6 * n1 + 2 * i + 1];
            }
        }
    } else if (n2 == 4 * nthread) {
        for (i = 0; i < n1; i++) {
            t[2 * i] = a[i][4 * n0];
            t[2 * i + 1] = a[i][4 * n0 + 1];
            t[2 * n1 + 2 * i] = a[i][4 * n0 + 2];
            t[2 * n1 + 2 * i + 1] = a[i][4 * n0 + 3];
        }
        cdft(2 * n1, isgn, t, ip, w);
        cdft(2 * n1, isgn, &t[2 * n1], ip, w);
        for (i = 0; i < n1; i++) {
            a[i][4 * n0] = t[2 * i];
            a[i][4 * n0 + 1] = t[2 * i + 1];
            a[i][4 * n0 + 2] = t[2 * n1 + 2 * i];
            a[i][4 * n0 + 3] = t[2 * n1 + 2 * i + 1];
        }
    } else if (n2 == 2 * nthread) {
        for (i = 0; i < n1; i++) {
            t[2 * i] = a[i][2 * n0];
            t[2 * i + 1] = a[i][2 * n0 + 1];
        }
        cdft(2 * n1, isgn, t, ip, w);
        for (i = 0; i < n1; i++) {
            a[i][2 * n0] = t[2 * i];
            a[i][2 * n0 + 1] = t[2 * i + 1];
        }
    }
    return (void *) 0;
}


void *ddxt2d0_th(void *p)
{
    void ddct(int n, int isgn, double *a, int *ip, double *w);
    void ddst(int n, int isgn, double *a, int *ip, double *w);
    int nthread, n0, n1, n2, ics, isgn, *ip, i;
    double **a, *w;
    
    nthread = ((fft2d_arg_t *) p)->nthread;
    n0 = ((fft2d_arg_t *) p)->n0;
    n1 = ((fft2d_arg_t *) p)->n1;
    n2 = ((fft2d_arg_t *) p)->n2;
    ics = ((fft2d_arg_t *) p)->ic;
    isgn = ((fft2d_arg_t *) p)->isgn;
    a = ((fft2d_arg_t *) p)->a;
    ip = ((fft2d_arg_t *) p)->ip;
    w = ((fft2d_arg_t *) p)->w;
    if (ics == 0) {
        for (i = n0; i < n1; i += nthread) {
            ddct(n2, isgn, a[i], ip, w);
        }
    } else {
        for (i = n0; i < n1; i += nthread) {
            ddst(n2, isgn, a[i], ip, w);
        }
    }
    return (void *) 0;
}


void *ddxt2d_th(void *p)
{
    void ddct(int n, int isgn, double *a, int *ip, double *w);
    void ddst(int n, int isgn, double *a, int *ip, double *w);
    int nthread, n0, n1, n2, ics, isgn, *ip, i, j;
    double **a, *t, *w;
    
    nthread = ((fft2d_arg_t *) p)->nthread;
    n0 = ((fft2d_arg_t *) p)->n0;
    n1 = ((fft2d_arg_t *) p)->n1;
    n2 = ((fft2d_arg_t *) p)->n2;
    ics = ((fft2d_arg_t *) p)->ic;
    isgn = ((fft2d_arg_t *) p)->isgn;
    a = ((fft2d_arg_t *) p)->a;
    t = ((fft2d_arg_t *) p)->t;
    ip = ((fft2d_arg_t *) p)->ip;
    w = ((fft2d_arg_t *) p)->w;
    if (n2 > 2 * nthread) {
        for (j = 4 * n0; j < n2; j += 4 * nthread) {
            for (i = 0; i < n1; i++) {
                t[i] = a[i][j];
                t[n1 + i] = a[i][j + 1];
                t[2 * n1 + i] = a[i][j + 2];
                t[3 * n1 + i] = a[i][j + 3];
            }
            if (ics == 0) {
                ddct(n1, isgn, t, ip, w);
                ddct(n1, isgn, &t[n1], ip, w);
                ddct(n1, isgn, &t[2 * n1], ip, w);
                ddct(n1, isgn, &t[3 * n1], ip, w);
            } else {
                ddst(n1, isgn, t, ip, w);
                ddst(n1, isgn, &t[n1], ip, w);
                ddst(n1, isgn, &t[2 * n1], ip, w);
                ddst(n1, isgn, &t[3 * n1], ip, w);
            }
            for (i = 0; i < n1; i++) {
                a[i][j] = t[i];
                a[i][j + 1] = t[n1 + i];
                a[i][j + 2] = t[2 * n1 + i];
                a[i][j + 3] = t[3 * n1 + i];
            }
        }
    } else if (n2 == 2 * nthread) {
        for (i = 0; i < n1; i++) {
            t[i] = a[i][2 * n0];
            t[n1 + i] = a[i][2 * n0 + 1];
        }
        if (ics == 0) {
            ddct(n1, isgn, t, ip, w);
            ddct(n1, isgn, &t[n1], ip, w);
        } else {
            ddst(n1, isgn, t, ip, w);
            ddst(n1, isgn, &t[n1], ip, w);
        }
        for (i = 0; i < n1; i++) {
            a[i][2 * n0] = t[i];
            a[i][2 * n0 + 1] = t[n1 + i];
        }
    } else if (n2 == nthread) {
        for (i = 0; i < n1; i++) {
            t[i] = a[i][n0];
        }
        if (ics == 0) {
            ddct(n1, isgn, t, ip, w);
        } else {
            ddst(n1, isgn, t, ip, w);
        }
        for (i = 0; i < n1; i++) {
            a[i][n0] = t[i];
        }
    }
    return (void *) 0;
}
#endif /* USE_FFT2D_THREADS */


//[end]

#endif