static void ccall(JF, js_Ast *fun, js_Ast *args)
{
	int n;
	switch (fun->type) {
	case EXP_INDEX:
		cexp(J, F, fun->a);
		emit(J, F, OP_DUP);
		cexp(J, F, fun->b);
		emit(J, F, OP_GETPROP);
		emit(J, F, OP_ROT2);
		break;
	case EXP_MEMBER:
		cexp(J, F, fun->a);
		emit(J, F, OP_DUP);
		emitstring(J, F, OP_GETPROP_S, fun->b->string);
		emit(J, F, OP_ROT2);
		break;
	case EXP_IDENTIFIER:
		if (!strcmp(fun->string, "eval")) {
			ceval(J, F, fun, args);
			return;
		}
		/* fall through */
	default:
		cexp(J, F, fun);
		emit(J, F, J->strict ? OP_UNDEF : OP_GLOBAL);
		break;
	}
	n = cargs(J, F, args);
	emit(J, F, OP_CALL);
	emitraw(J, F, n);
}

/* Conformal map that takes the unit disc onto the disc minus a slit
 * starting at cexp(I * u). Approximates the atomic maps for radial
 * SLE, and the reciprocal of this approximates the atomic maps for
 * full plane SLE. */
static complex sle_cmap_slit_disc(double t, double u, complex z)
{
  complex rotation = cexp((I * u));
  complex w = 4 * cexp(-t) * z;
  complex v = z + 1;
  return rotation * 1.0/w * (2 * v * v - w - 2 * v * csqrt(v * v - w));
}

void fft(double _Complex *in, double _Complex *out,unsigned char N, short sign)
{
	double _Complex *mid;
	unsigned temp, temp1,rank;
	unsigned listlenth = pow(2,N);

	mid=(double _Complex *)malloc(sizeof(double _Complex)*listlenth);
	for(unsigned i=0;i<listlenth;i++)
		mid[i]=in[reversebit(i,(unsigned char)N)];

	unsigned j, t, interval;
	double _Complex exponent, exponent_unit;

	interval = 1;
	for(unsigned i=1;i<=N;i++){
		j=0;
		while(j<listlenth){
			exponent_unit = (sign<=0)? cexp(I*PI/interval):cexp(-I*PI/interval);
			exponent= 1;
			for(t=0;t<interval*2;t++){
				out[t+j]=mid[t%interval+j]+exponent*mid[t%interval+interval+ j];
				exponent *= exponent_unit;
			}
			j=j+interval*2;
		}
		interval *= 2;
		memcpy(mid,out,listlenth*sizeof(double _Complex));
	}

	if (sign <=0){
		for(j = 0;j <listlenth;j++)
			out[j] /= listlenth;
	}
	free(mid);
}

LTFAT_EXTERN void
LTFAT_NAME_REAL(idwiltiii_long)(const LTFAT_REAL *c, const LTFAT_REAL *g,
                                const ltfatInt L, const ltfatInt W,
                                const ltfatInt M, LTFAT_REAL *f)
{
    const ltfatInt N = L / M;
    const ltfatInt M2 = 2 * M;
    const ltfatInt M4 = 4 * M;
    const LTFAT_REAL scalconst = 1.0 / sqrt(2.0);
    const LTFAT_COMPLEX eipi4 = cexp(I * PI / 4.0);
    const LTFAT_COMPLEX emipi4 = cexp(-I * PI / 4.0);

    LTFAT_COMPLEX *coef2 = ltfat_calloc(2 * M * N * W, sizeof * coef2);
    LTFAT_COMPLEX *f2 = ltfat_malloc(L * W * sizeof * f2);
    LTFAT_COMPLEX *g2 = ltfat_malloc(L * sizeof * g2);
    for (ltfatInt ii = 0; ii < L; ii++)
        g2[ii] = g[ii];


    const LTFAT_REAL *pcoef  = c;
    LTFAT_COMPLEX *pcoef2 = coef2;

    PREPROC_COMPLEX

    LTFAT_NAME(idgt_long)(coef2, g2, L, W, M, 2 * M, FREQINV, f2);

    POSTPROC_REAL

    LTFAT_SAFEFREEALL(coef2, f2, g2);

}

LTFAT_EXTERN void
LTFAT_NAME_COMPLEX(dwiltiii_long)(const LTFAT_COMPLEX *f, const LTFAT_COMPLEX *g,
                                  const ltfatInt L, const ltfatInt W, const ltfatInt M,
                                  LTFAT_COMPLEX *cout)
{
    const ltfatInt N=L/M;
    const ltfatInt M2=2*M;
    const ltfatInt M4=4*M;
    const LTFAT_REAL scalconst = 1.0/sqrt(2.0);
    const LTFAT_COMPLEX eipi4 = cexp(I*PI/4.0);
    const LTFAT_COMPLEX emipi4 = cexp(-I*PI/4.0);

    LTFAT_COMPLEX *coef2 = ltfat_malloc(2*M*N*W*sizeof*coef2);
    LTFAT_COMPLEX *f2 = ltfat_malloc(L*W*sizeof*f2);

    PREPROC

    LTFAT_NAME_COMPLEX(dgt_long)(f2, g, L, W, M, 2*M, coef2);

    LTFAT_COMPLEX *pcoef  = cout;
    LTFAT_COMPLEX *pcoef2 = coef2;

    POSTPROC_COMPLEX

    LTFAT_SAFEFREEALL(coef2,f2);
}

static void ccall(JF, js_Ast *fun, js_Ast *args)
{
	int n;
	switch (fun->type) {
	case EXP_INDEX:
		cexp(J, F, fun->a);
		emit(J, F, OP_DUP);
		cexp(J, F, fun->b);
		emit(J, F, OP_GETPROP);
		emit(J, F, OP_ROT2);
		break;
	case EXP_MEMBER:
		cexp(J, F, fun->a);
		emit(J, F, OP_DUP);
		emitstring(J, F, OP_GETPROP_S, fun->b->string);
		emit(J, F, OP_ROT2);
		break;
	default:
		cexp(J, F, fun);
		emit(J, F, OP_GLOBAL);
		break;
	}
	n = cargs(J, F, args);
	emit(J, F, OP_CALL);
	emitraw(J, F, n);
}

void mkexps(const double *rlams, const int nlambs, const int *numphys,
	    const int nexptotp, dcomplex *xs, dcomplex *ys, double *zs)
{
  int ntot, nell, mth, ncurrent;
  double hu, u;

  ntot = 0;
  for ( nell = 0; nell < nlambs; nell++) {
    hu = 2.0*M_PI/numphys[nell];
    for ( mth = 0; mth < numphys[nell]/2; mth++ ) {
      u = mth*hu;
      ncurrent=3*(ntot+mth);
      zs[ncurrent] = exp( -rlams[nell] );
      zs[ncurrent+1] = zs[ncurrent]*zs[ncurrent];
      zs[ncurrent+2] = zs[ncurrent]*zs[ncurrent+1];
      xs[ncurrent] = cexp(_Complex_I*cos(u)*rlams[nell]);
      xs[ncurrent+1] = xs[ncurrent]*xs[ncurrent];
      xs[ncurrent+2] = xs[ncurrent]*xs[ncurrent+1];
      ys[ncurrent] = cexp(_Complex_I*sin(u)*rlams[nell]);
      ys[ncurrent+1] = ys[ncurrent]*ys[ncurrent];
      ys[ncurrent+2] = ys[ncurrent+1]*ys[ncurrent];
    }
    ntot += numphys[nell]/2;
  }
}

inline  Complex PlasmaDispersion(Complex z) {
   
        static const double  SQRT_PI = 1.7724538509055160272981674833411451;
        Complex  Z;
        Complex z2 = z*z;
       

        // sometimes imaginary part is too high, due to test functions, limit
        if(cimag(z) > 25.) z = creal(z) +  25.I;

        // Use asymptotic expansions for large values |z| >> 1
        // See definition in NRL Plasma Formulary
        if(creal(z*z) > 100) {
            
            const double x = creal(z); 
            const double y = cimag(z);
            
            const double sigma = (y >  0.) ? 0 : (( y == 0) ? 1. : 2.);
            Complex z4 = z2*z2;
            Complex z6 = z4*z2;
            Z = 1.0I *   SQRT_PI * sigma * cexp(-z2) - 1./z * ( 1. + 1./(2.*z2) + 3./(4.*z4) + 15./(8.*z6));		
  
        }
        else Z = 1.0I * SQRT_PI * cexp(- z2) * cerfc(- 1.0I * z);
    return Z;
}; 

LTFAT_EXTERN void
LTFAT_NAME_COMPLEX(idwiltiii_fb)(const LTFAT_COMPLEX *c, const LTFAT_COMPLEX *g,
                                 const ltfatInt L, const ltfatInt gl,
                                 const ltfatInt W, const ltfatInt M,
                                 LTFAT_COMPLEX *f)
{
    const ltfatInt N = L / M;
    const ltfatInt M2 = 2 * M;
    const ltfatInt M4 = 4 * M;
    const LTFAT_REAL scalconst = 1.0 / sqrt(2.0);
    const LTFAT_COMPLEX eipi4 = cexp(I * PI / 4.0);
    const LTFAT_COMPLEX emipi4 = cexp(-I * PI / 4.0);

    LTFAT_COMPLEX *coef2 = ltfat_calloc(2 * M * N * W, sizeof * coef2);
    LTFAT_COMPLEX *f2 = ltfat_malloc(L * W * sizeof * f2);


    const LTFAT_COMPLEX *pcoef  = c;
    LTFAT_COMPLEX *pcoef2 = coef2;

    PREPROC_COMPLEX

    LTFAT_NAME(idgt_fb)(coef2, g, L, gl, W, M, 2 * M, FREQINV, f2);

    POSTPROC_COMPLEX

    LTFAT_SAFEFREEALL(coef2, f2);

}

static void cassignforin(JF, js_Ast *stm)
{
	js_Ast *lhs = stm->a;

	if (stm->type == STM_FOR_IN_VAR) {
		if (lhs->b)
			jsC_error(J, lhs->b, "more than one loop variable in for-in statement");
		emitlocal(J, F, OP_SETLOCAL, OP_SETVAR, lhs->a->a); /* list(var-init(ident)) */
		emit(J, F, OP_POP);
		return;
	}

	switch (lhs->type) {
	case EXP_IDENTIFIER:
		emitlocal(J, F, OP_SETLOCAL, OP_SETVAR, lhs);
		emit(J, F, OP_POP);
		break;
	case EXP_INDEX:
		cexp(J, F, lhs->a);
		cexp(J, F, lhs->b);
		emit(J, F, OP_ROT3);
		emit(J, F, OP_SETPROP);
		emit(J, F, OP_POP);
		break;
	case EXP_MEMBER:
		cexp(J, F, lhs->a);
		emit(J, F, OP_ROT2);
		emitstring(J, F, OP_SETPROP_S, lhs->b->string);
		emit(J, F, OP_POP);
		break;
	default:
		jsC_error(J, lhs, "invalid l-value in for-in loop assignment");
	}
}

static void calculate(catastrophe_t *const catastrophe,
		const unsigned  int i, const unsigned int j)
{
	cmplx_equation_t *equation;
	point_array_t *point_array;

	/* Gamma function's precalculated values */
	const double g14 = 3.625609908;
	const double g34 = 1.225416702;

	assert(catastrophe);

	equation = catastrophe->equation;
	point_array = catastrophe->point_array;

	assert(equation);
	assert(point_array);

	equation->initial_vector[V] = 0.5 * g14 * cexp(I * M_PI / 8.0);
	equation->initial_vector[V1] = 0;
	equation->initial_vector[V2] = 0.5 * I * g34 *
		cexp(I * 3.0 * M_PI / 8.0);

	cmplx_runge_kutta(0.0, 1.0, 0.01, catastrophe);
}

/**
   Make basic arrays for servo analysis.
*/
static void servo_calc_init(SERVO_CALC_T *st, const dmat *psdin, double dt, long dtrat){
    if(psdin->ny!=2){
	error("psdin should have two columns\n");
    }
    double Ts=dt*dtrat;
    st->fny=0.5/Ts;
    dmat *nu=st->nu=dlogspace(-3,log10(0.5/dt),1000);
    st->psd=dinterp1(psdin, 0, nu, 1e-40);
    st->var_sig=psd_inte2(psdin);
    dcomplex pi2i=COMPLEX(0, TWOPI);
    if(st->Hsys || st->Hwfs || st->Hint || st->s){
	error("Already initialized\n");
    }
    st->Hsys=cnew(nu->nx, 1);
    st->Hwfs=cnew(nu->nx, 1);
    st->Hint=cnew(nu->nx, 1);
    st->s=cnew(nu->nx,1);
    for(long i=0; i<nu->nx; i++){
	dcomplex s=st->s->p[i]=pi2i*nu->p[i];
	dcomplex zInv=cexp(-s*Ts);
	dcomplex Hint=st->Hint->p[i]=1./(1-zInv);
	dcomplex Hwfs, Hdac;
	if(dtrat==1){//we have a pure delay
	    Hwfs=st->Hwfs->p[i]=zInv; 
	    Hdac=1;
	}else{
	    Hwfs=st->Hwfs->p[i]=(1-zInv)/(Ts*s);
	    Hdac=Hwfs;
	}
	dcomplex Hlag=cexp(-s*dt);/*lag due to readout/computation*/
	dcomplex Hmir=1;/*DM */
	st->Hsys->p[i]=Hwfs*Hlag*Hdac*Hmir*Hint;
    }
}

/**
 * Initialisation of a transform plan, guru.
 *
 * \arg ths The pointer to a fpt plan
 * \arg N The number of source nodes
 * \arg M The number of target nodes
 * \arg sigma The parameter of the Gaussian
 * \arg n The polynomial expansion degree
 * \arg p the periodisation length, at least 1
 * \arg m The spatial cut-off of the nfft
 * \arg flags FGT flags to use
 *
 * \author Stefan Kunis
 */
void fgt_init_guru(fgt_plan *ths, int N, int M, double _Complex sigma, int n,
		   double p, int m, unsigned flags)
{
  int j,n_fftw;
  fftw_plan fplan;

  ths->M = M;
  ths->N = N;
  ths->sigma = sigma;
  ths->flags = flags;

  ths->x = (double*)nfft_malloc(ths->N*sizeof(double));
  ths->y = (double*)nfft_malloc(ths->M*sizeof(double));
  ths->alpha = (double _Complex*)nfft_malloc(ths->N*sizeof(double _Complex));
  ths->f = (double _Complex*)nfft_malloc(ths->M*sizeof(double _Complex));

  ths->n = n;
  ths->p = p;

  ths->b = (double _Complex*)nfft_malloc(ths->n*sizeof(double _Complex));

  ths->nplan1 = (nfft_plan*) nfft_malloc(sizeof(nfft_plan));
  ths->nplan2 = (nfft_plan*) nfft_malloc(sizeof(nfft_plan));

  n_fftw=X(next_power_of_2)(2*ths->n);

  nfft_init_guru(ths->nplan1, 1, &(ths->n), ths->N, &n_fftw, m, PRE_PHI_HUT|
                 PRE_PSI| MALLOC_X| MALLOC_F_HAT| FFTW_INIT, FFTW_MEASURE);
  nfft_init_guru(ths->nplan2, 1, &(ths->n), ths->M, &n_fftw, m, PRE_PHI_HUT|
                 PRE_PSI| MALLOC_X| FFTW_INIT, FFTW_MEASURE);

  ths->nplan1->f = ths->alpha;
  ths->nplan2->f_hat = ths->nplan1->f_hat;
  ths->nplan2->f = ths->f;

  if(ths->flags & FGT_APPROX_B)
    {
      fplan = fftw_plan_dft_1d(ths->n, ths->b, ths->b, FFTW_FORWARD,
                               FFTW_MEASURE);

      for(j=0; j<ths->n; j++)
	ths->b[j] = cexp(-ths->p*ths->p*ths->sigma*(j-ths->n/2)*(j-ths->n/2)/
                          ((double)ths->n*ths->n)) / ths->n;

      nfft_fftshift_complex(ths->b, 1, &ths->n);
      fftw_execute(fplan);
      nfft_fftshift_complex(ths->b, 1, &ths->n);

      fftw_destroy_plan(fplan);
    }
  else
    {
      for(j=0; j<ths->n; j++)
	ths->b[j] = 1.0/ths->p * csqrt(PI/ths->sigma)*
	  cexp(-PI*PI*(j-ths->n/2)*(j-ths->n/2)/
	       (ths->p*ths->p*ths->sigma));
    }
}

/*! \memberof splitop
 create new operator using values in prefs */
splitop_t * splitop_new(preferences_t * prefs)
{
  splitop_t * w = malloc(sizeof(splitop_t)); assert(w);
  w->prefs = prefs;

  // cache values
  double * V = prefs->potential->data;
  int bins = prefs->bins;
  double dt = prefs->dt;
  double dk = prefs->dk;

  // allocate all the arrays needed
  w->eV = fftw_alloc_complex(bins); assert(w->eV);
  w->eVn = fftw_alloc_complex(bins); assert(w->eVn);
  w->ehV = fftw_alloc_complex(bins); assert(w->ehV);
  w->ehVn = fftw_alloc_complex(bins); assert(w->ehVn);
  w->eT = fftw_alloc_complex(bins); assert(w->eT);
  w->apsi = fftw_alloc_complex(bins); assert(w->apsi);
  w->psi = fftw_alloc_complex(bins); assert(w->psi);
  w->psik = fftw_alloc_complex(bins); assert(w->psik);

  for (int k = 0; k < bins; k++) {
    w->apsi[k] = 0;
  }

  // calculate position space propagators
  for (int n = 0; n < bins; n++) {
    w->eV[n] = cexp(-I * V[n] * dt);
    w->eVn[n] = w->eV[n] / bins;
    w->ehV[n] = cexp(-I * V[n] * dt / 2);
    w->ehVn[n] = w->ehV[n] / bins;
  }

  // calculate momentum space propagator
  for (int n = 0; n < bins; n++) {
    if (n < bins/2) {
      w->eT[n] = cexp(-I * dk * dk * n * n * dt);
    } else { // negative frequencies come after bins/2
      int m = (n - bins);
      w->eT[n] = cexp(-I * dk * dk * m * m * dt);
    }
  }


  // prepare fftw plans for DFT
  w->fwd = fftw_plan_dft_1d(bins, w->psi, w->psik, FFTW_FORWARD, FFTW_MEASURE);
  w->bwd = fftw_plan_dft_1d(bins, w->psik, w->psi, FFTW_BACKWARD, FFTW_MEASURE);

  // setup initial system state wavefunction
  for (int k = 0; k < bins; k++) {
    w->psi[k] = prefs->psi->data[k];
  }

  return w;
}

static complex     phi_recov_unif(const double     u, const void         *p_recovery)
{
  params_recov_unif   *p = (params_recov_unif *) p_recovery;
  complex     z1;
  complex     z2;
  complex         result;

  if (u == 0) return(1.);
  z1 = cexp(I*u*(1-p->a)) - cexp(I*u*(1-p->b));
  z2 = I*u*(p->b-p->a);
  result = z1 / z2;
  return ( result );
}

void init_LargeScale2DNoise(struct Field fldi) {
	int i,j,k;
	int num_force=0;
	int total_num_force;
	double fact;
	
	for( i = 0; i < NX_COMPLEX/NPROC; i++) {
		for( j = 0; j < NY_COMPLEX; j++) {
			k=0;
			if(kz[ IDX3D ] == 0.0) {
				if(pow(k2t[ IDX3D ], 0.5) / ( 2.0*M_PI ) < 1.0 / param.noise_cut_length_2D) {
					fldi.vx[ IDX3D ] += param.per_amplitude_large_2D * mask[IDX3D] * randm() * cexp( I * 2.0*M_PI*randm() ) * NTOTAL;
					fldi.vy[ IDX3D ] += param.per_amplitude_large_2D * mask[IDX3D] * randm() * cexp( I * 2.0*M_PI*randm() ) * NTOTAL;
#ifdef MHD
					fldi.bx[ IDX3D ] += param.per_amplitude_large_2D * mask[IDX3D] * randm() * cexp( I * 2.0*M_PI*randm() ) * NTOTAL;
					fldi.by[ IDX3D ] += param.per_amplitude_large_2D * mask[IDX3D] * randm() * cexp( I * 2.0*M_PI*randm() ) * NTOTAL;
#endif
					if(mask[IDX3D] > 0) num_force++;
				}
			}
		}
	}
	
	// Get the total number of forced scales.
#ifdef MPI_SUPPORT
	MPI_Allreduce( &num_force, &total_num_force, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
#else
	total_num_force=num_force;
#endif
	
	fact=pow(total_num_force,0.5);
	
	// Divide by the total number of modes
	for( i = 0; i < NX_COMPLEX/NPROC; i++) {
		for( j = 0; j < NY_COMPLEX; j++) {
			k=0;
			if(kz[ IDX3D ] == 0.0) {
				fldi.vx[ IDX3D ] = fldi.vx[ IDX3D ] / fact;
				fldi.vy[ IDX3D ] = fldi.vy[ IDX3D ] / fact;
#ifdef MHD
				fldi.bx[ IDX3D ] = fldi.bx[ IDX3D ] / fact;
				fldi.by[ IDX3D ] = fldi.by[ IDX3D ] / fact;
#endif
			}
		}
	}
	
  enforce_complex_symm(fldi);  
}

static void findroots(complex double a, complex double b, complex double c,
							complex double e[])
{
	int i;
	complex double p = (b - a * a / 3.) / 3.;
	complex double q = (2. * a * a * a / 27. - a * b / 3. + c) / 2.;
	complex double D = creal(p * p * p + q * q); /* MODEL SPECIFIC!!! */
	double eps;

	assert(abs(cimag(p)) < 1e-10);
	assert(abs(creal(q)) < 1e-10);
	assert(abs(cimag(p * p * p + q * q)) < 1e-10);

	complex double upr = -q + csqrt(D);
	double mod_upr = pow(cabs(upr), 1. / 3.);
	double arg_upr = carg(upr);
	complex double umr = -q - csqrt(D);
	double mod_umr = pow(cabs(umr), 1. / 3.);
	double arg_umr = carg(umr);

	complex double rp = .5 * (-1. + I * sqrt(3.));
	complex double rm = .5 * (-1. - I * sqrt(3.));
	complex double up = mod_upr * cexp(I * arg_upr / 3.);

	for (eps = 1e-30; eps < 1e-6; eps *= 10.) {
		complex double um[3];
		double sort[3];

		for (i = 0; i < 3; i++) {
			um[i] = mod_umr * cexp(I*(2. * M_PI * i + arg_umr) / 3.);
			sort[i] = cabs((um[i] * up + p) / p);
		}
		qsort(sort, 3, sizeof(*sort), cmp);

		for (i = 0; i < 3; i++) {
			double test = cabs((um[i] * up + p) / p);
			if (test == sort[0] && test < eps) {
				e[0] = up + um[i] - a / 3.;
				e[1] = rp * up + rm * um[i] - a / 3.;
				e[2] = rm * up + rp * um[i] - a / 3.;
				return;
			}
		}
	}
	fprintf(stderr, "This should never happen!\n");
	fprintf(stderr, "up=%lg%+lg*I p=%lg%+lg*I q=%lg%+lg*I D=%lg\n",
		creal(up), cimag(up), creal(p), cimag(p),
		creal(q), cimag(q), creal(D));
}

//! Creates the test signal
// 0 (default): random
// 1 : lin. comb. of exp's
// 2 : exp + noise
// 3 : sinc
// 4 : gaussian
// 5 : audio
void create_signal(complex double *vec, int sz, int kind){
	int i, sz2;
	double tmp;
	double *buf;
	FILE *fin;

	sz2 = sz>>1;
	switch (kind){
		case 1:
			// sum of exp's -- tiled version
			for (i=0; i<sz; i++) vec[i] = cexp(I*i*M_PI*0.05) + exp(I*i*M_PI*0.075);
			break;
		case 2:
			// exp + noise -- tiled version
			srand(0);
			for (i=0; i<sz; i++) vec[i] = cexp(I*i*M_PI*0.05) + (rand()%(2001) - 1000)*1e-4;
			break;
		case 3:
			// sinc -- rescaled version
			tmp = 100./(double) sz;
			for (i=0; i<sz; i++) vec[i] = sin(tmp*(i-sz+1e-4)) / (tmp*(i-sz+1e-4));
			break;
		case 4:
			// gaussian -- rescalled version
			for (i=0; i<sz; i++) vec[i] = exp(-tmp*pow(i-sz2,2));
			break;
		case 5:
			// audio -- fixed length
			assert(sizeof(double)==8);
			fin = fopen("tests/music.dat", "rb");
			if (!fin) {
				printf("Can't open tests/music.dat. Make sure the file is there!\n");
				return;
			}
			buf = malloc( (1<<17)*sizeof(double) );
			//jump initial noisy part
			fread(buf, sizeof(double), (1<<12) , fin);
			fread(buf, sizeof(double), (1<<17) , fin);
			fclose(fin);
			for (i=0; i<sz; i++) vec[i] = (double complex) buf[i%(1<<17)];
			free(buf);
			break;
		default:
			srand(0);
			for (i=0; i<sz; i++) vec[i] = rand()%2001 - 1000;
			break;
	}
}

Gamma * gamma_spline_calc_gamma (const GammaSpline *gs, double epsilon) {
	//int q=0;
	double kc;
    int Nkc = Nkc_from_epsilon(epsilon);
    Gamma * g = gamma_new (Nkc);
    double dkc=2*KCMAX/Nkc;
    double denom=pow(KCMAX,4);
    
    int q=Nkc/2;   /* we do this in a funny order to keep phase oscillations out of the spectrum */
	for (kc=-KCMAX+dkc;kc<KCMAX;kc=kc+dkc) {
	    complex double a = gsl_interp_eval(gs->wx_spline_re,gs->k,gs->Wxr,kc,gs->w_accel)+I*gsl_interp_eval(gs->wx_spline_im,gs->k,gs->Wxi,kc,gs->w_accel);
	    complex double b = cexp(-4*pow(fabs(kc),4)/denom+I*gsl_interp_eval(gs->phi_spline,gs->k,gs->phi,kc,gs->phi_accel)/epsilon);
	    g->x[q]= a*b;
	    
	    a = gsl_interp_eval(gs->wy_spline_re,gs->k,gs->Wyr,kc,gs->w_accel)+I*gsl_interp_eval(gs->wy_spline_im,gs->k,gs->Wyi,kc,gs->w_accel);
	    g->y[q]=a*b;
	    
	    a = gsl_interp_eval(gs->wz_spline_re,gs->k,gs->Wzr,kc,gs->w_accel)+I*gsl_interp_eval(gs->wz_spline_im,gs->k,gs->Wzi,kc,gs->w_accel);
	    
	    g->z[q]=a*b;
	    
	    q++;
	    if (q==Nkc)
	        q=0;
	}
	return g;
}

// Function to generate the DFT/IDFT N x N Twiddle matrix to be passed to the GPU
//! @param isIDFT 				IDFT/DFT option
//! @param pTwiddleMatrix		Generated twiddle matrix handle
//! @param pTwiddleMatrix_z		Generated twiddle matrix complex type handle
void twiddleMatrixGen(bool isIDFT, Complex *pTwiddleMatrix, double complex *pTwiddleMatrix_z)
{
	// Exponent return value
	double complex m_result;

	// Capture Twiddle Matrix reference for resetting pointer before return
	Complex *pTwiddleRef = pTwiddleMatrix;
	double complex *pTwiddleRef_z = pTwiddleMatrix_z;

	// Condition exponent sign bit and scaling on DFT or IDFT
	int sign = -1;
	float scale = 1;
	if ( !isIDFT )
	{
		sign = 1;
		scale = 1 / static_cast<float> ( N ) ;
	}

	for ( int i = 0; i < N; ++i )
		for ( int j = 0; j < N; ++j, ++pTwiddleMatrix, ++pTwiddleMatrix_z )
		{
			m_result = cexp( sign * 2 * M_PI * i * j / N * I ) * scale;
			*pTwiddleMatrix_z = m_result;
			pTwiddleMatrix->x = creal( m_result );
			pTwiddleMatrix->y = cimag( m_result );
		}

	// Reset the Matrix pointers
	pTwiddleMatrix = pTwiddleRef;
	pTwiddleMatrix_z = pTwiddleRef_z;
}