main()
{
	int npoles,nfft,i;
	float f3db,zero=0.0,fpass,apass,fstop,astop;

/*
	printf("Enter npoles f3db:\n");
	scanf("%d %f",&npoles,&f3db);
*/
	printf("Enter fpass apass fstop astop:\n");
	scanf("%f %f %f %f",&fpass,&apass,&fstop,&astop);
	bfdesign(fpass,apass,fstop,astop,&npoles,&f3db);
	printf("npoles = %d  f3db = %f\n",npoles,f3db);

	/* impulse response */
	scopy(N,&zero,0,p,1);
	p[0] = 1.0;
	bflowpass(npoles,f3db,N,p,q);
	pp1d(stdout,"impulse response",N,0,q);

	/* amplitude spectrum */
	nfft = npfa(N);
	for (i=0; i<N; i++)
		z[i] = cmplx(q[i],0.0);
	for (i=N; i<nfft; i++)
		z[i] = cmplx(0.0,0.0);
	pfacc(1,nfft,z);
	for (i=0; i<nfft; i++)
		zamp[i] = fcabs(z[i]);
	pp1d(stdout,"amplitude spectrum",nfft/2+1,0,zamp);
}

void CIIRfilter::butter_poles(complex *p, char *type, int *n, int iord )
	{
	double	pi, angle;
	int	half, k;

        pi=M_PI;
        half = iord/2;

	/* test for odd order, and add pole at -1 */

        *n = 0;
        if (    2*half < iord )
			{
			p[0] = cmplx( -1., 0. );
			type[0] = 'S';
			*n = 1;
			}
	for (k = 0; k < half; k++)
		{
		angle = pi * ( 0.5f + (double)(2*(k+1)-1) / (double)(2*iord) );
		p[*n] = cmplx( (double)cos((double)angle), 
			(double)sin((double)angle) );
		type[*n] = 'C';
		*n = *n + 1;
		}
}

/* do measurements: load density, ploop, etc. and phases onto lattice */
void measure() {
   register int i,j,k, c, is_even;
   register site *s;
   int dx,dy,dz;	/* separation for correlated observables */
   int dir;		/* direction of separation */
   msg_tag *tag;
   register complex cc,dd;	/*scratch*/
   complex ztr, zcof, znum, zdet, TC, zd, density, zphase;
   complex p[4]; /* probabilities of n quarks at a site */
   complex np[4]; /* probabilities at neighbor site */
   complex pp[4][4]; /* joint probabilities of n here and m there */
   complex zplp, plp_even, plp_odd;
   Real locphase, phase;


   /* First make T (= timelike P-loop) from s->ploop_t 
      T stored in s->tempmat1
   */
   ploop_less_slice(nt-1,EVEN);
   ploop_less_slice(nt-1,ODD);

   phase = 0.;
   density = plp_even = plp_odd = cmplx(0.0, 0.0);
   for(j=0;j<4;j++){
	p[j]=cmplx(0.0,0.0);
	for(k=0;k<4;k++)pp[j][k]=cmplx(0.0,0.0);
   }
   FORALLSITES(i,s) {
      if(s->t != nt-1) continue;
      if( ((s->x+s->y+s->z)&0x1)==0 ) is_even=1; else is_even=0;
      mult_su3_nn(&(s->link[TUP]), &(s->ploop_t), &(s->tempmat1));

      zplp = trace_su3(&(s->tempmat1));
      if( is_even){CSUM(plp_even, zplp)}
      else        {CSUM(plp_odd, zplp)}

      ztr = trace_su3(&(s->tempmat1));
      CONJG(ztr, zcof);

      if(is_even){
        for(c=0; c<3; ++c) s->tempmat1.e[c][c].real += C;
        zdet = det_su3(&(s->tempmat1));
        znum = numer(C, ztr, zcof);
        CDIV(znum, zdet, zd);
        CSUM(density, zd);

        /* store n_quark probabilities at this site in lattice variable
	  qprob[], accumulate sum over lattice in p[] */
        cc= cmplx(C*C*C,0.0); CDIV(cc,zdet,s->qprob[0]); CSUM(p[0],s->qprob[0]);
        CMULREAL(ztr,C*C,cc); CDIV(cc,zdet,s->qprob[1]); CSUM(p[1],s->qprob[1]);
        CMULREAL(zcof,C,cc); CDIV(cc,zdet,s->qprob[2]); CSUM(p[2],s->qprob[2]);
        cc = cmplx(1.0,0.0); CDIV(cc,zdet,s->qprob[3]); CSUM(p[3],s->qprob[3]);
  
        locphase = carg(&zdet);
        phase += locphase;
      }

   }

void relax(int NumStp)
{
  /* Do overrelaxation by SU(2) subgroups */
int NumTrj, Nhit, index1, ina, inb, ii;
int subl;
Real a0, a1, a2, a3, asq, r;
register int dir, i;
register site *st;
su3_matrix action;
su2_matrix u;

Nhit = 3;

    for( NumTrj = 0 ; NumTrj < NumStp; NumTrj++)
    for( subl = 0; subl < N_SUBL32; subl++) {

	FORALLUPDIR(dir) {
	    /* updating links in direction dir */
	    /* compute the staple */
	    dsdu_qhb_subl(dir, subl);

	    for(index1=0;index1< Nhit;index1++) {
		/*  pick out an SU(2) subgroup */
		ina=(index1+1) % Nc;
		inb=(index1+2) % Nc;
		if(ina > inb){ ii=ina; ina=inb; inb=ii;}

		FORSOMESUBLATTICE(i,st,subl) {
		    mult_su3_na(&(st->link[dir]), &(st->staple), &action);

/*decompose the action into SU(2) subgroups using Pauli matrix expansion */
/* The SU(2) hit matrix is represented as a0 + i * Sum j (sigma j * aj)*/
		    a0 =  action.e[ina][ina].real + action.e[inb][inb].real;
		    a3 =  action.e[ina][ina].imag - action.e[inb][inb].imag;
		    a1 =  action.e[ina][inb].imag + action.e[inb][ina].imag;
		    a2 =  action.e[ina][inb].real - action.e[inb][ina].real;

		    /* Normalize and complex conjugate u */

		    asq = a0*a0 + a1*a1 + a2*a2 + a3*a3;
		    r = sqrt((double)asq );
		    a0 = a0/r; a1 = -a1/r; a2 = -a2/r; a3 = -a3/r;

		    /* Elements of SU(2) matrix */
		    u.e[0][0] = cmplx( a0, a3);
		    u.e[0][1] = cmplx( a2, a1);
		    u.e[1][0] = cmplx(-a2, a1);
		    u.e[1][1] = cmplx( a0,-a3);

		    /* Do SU(2) hit on all links twice (to overrelax)  */
		    left_su2_hit_n(&u,ina,inb,&(st->link[dir]));
		    left_su2_hit_n(&u,ina,inb,&(st->link[dir]));

		} /*  st */
	    } /*   hits */
	} /*  direction */
    } /*  subl, NumTrj */

void CIIRfilter::lp_to_br(complex *p, char *ptype, int np, double fl, double fh, double *sn, double *sd, int *ns )
{
	complex	pinv, ctemp, p1, p2;
	double	pi, twopi, a, b;
	int	i, iptr;
        pi = M_PI;
        twopi = 2.*pi;
        a = twopi*twopi*fl*fh;
        b = twopi*( fh - fl );
        *ns = 0;
        iptr = 0;
	for (i = 0; i < np; i++)
	{
		if (    ptype[i] == 'C' )
		{
			pinv = cmplx_div(cmplx(1.,0.), p[i]);
			ctemp = cmplx_mul(real_cmplx_mul(b,pinv), 
				real_cmplx_mul(b, pinv));
			ctemp = cmplx_sub(ctemp, cmplx(4.*a, 0.));
			ctemp = cmplx_sqrt( ctemp );
			p1 = real_cmplx_mul(0.5, 
				cmplx_add(real_cmplx_mul(b,pinv), ctemp));
			p2 = real_cmplx_mul(0.5, 
				cmplx_sub(real_cmplx_mul(b,pinv), ctemp));
			sn[ iptr ]     = a;
			sn[ iptr + 1 ] = 0.;
			sn[ iptr + 2 ] = 1.;
			sd[iptr] = real_part(cmplx_mul(p1, cmplx_conjg(p1)));
			sd[ iptr + 1 ] = -2. * real_part(p1);
			sd[ iptr + 2 ] = 1.;
			iptr = iptr + 3;
			sn[ iptr ]     = a;
			sn[ iptr + 1 ] = 0.;
			sn[ iptr + 2 ] = 1.;
			sd[iptr] = real_part(cmplx_mul(p2, cmplx_conjg(p2)));
			sd[ iptr + 1 ] = -2. * real_part(p2);
			sd[ iptr + 2 ] = 1.;
			iptr = iptr + 3;
			*ns = *ns + 2;
		}
		else if (  ptype[i] == 'S' )
		{
			sn[ iptr ]     = a;
			sn[ iptr + 1 ] = 0.;
			sn[ iptr + 2 ] = 1.;
			sd[ iptr ]     = -a*real_part( p[i] );
			sd[ iptr + 1 ] = b;
			sd[ iptr + 2 ] = -real_part( p[i] );
			iptr = iptr + 3;
			*ns = *ns + 1;
		}
	}
}

/************************************************************************
dopow - raise a complex number to seperate real powers for amp and phase
*************************************************************************
Notes:

Aout exp[-i Pout] = Ain^a exp[-i b Pin] 
     where A=amplitude  P=phase
*************************************************************************
UHouston: Chris Liner: modified from D. Hale subroutine crpow
*************************************************************************/
complex dopow(complex u, float a, float b)
{
	float ur,ui,amp,phs;

	if (a==0.0) return cmplx(1.0,0.0);
	if (u.r==0.0 && u.i==0.0) return cmplx(0.0,0.0);

	ur = u.r; 
	ui = u.i;
	amp = exp(0.5*a*log(ur*ur+ui*ui));
	phs = -b*atan2(ui,ur);

	return cmplx(amp*cos(phs),amp*sin(phs));	
}

/******************* Wait for buffer of data to be input/output **********************/
void wait_buffer(void)
{
	float *p;  
	int i;
	/* wait for array index to be set to zero by ISR */
	while(index);
	
	/* rotate data arrays */
	p = input;
	input = output;
	output = intermediate;
	intermediate = p;
	
	/************************* DO PROCESSING OF FRAME  HERE **************************/                
	for(i=0;i<BUFLEN;i++){
		buffer[i]=cmplx(intermediate[i],0);
	}
	
	fft(BUFLEN,buffer);
	
	for(i=0;i<BUFLEN;i++){
		mag[i]=cabs(buffer[i]);
	}
	
	/**********************************************************************************/
	
	/* wait here in case next sample has not yet been read in */                          
	while(!index);
}        

void LapH::Correlators::compute_correlators(const size_t config_i){

  // initialising the big arrays
  set_corr(config_i);
  // setting the correlation functions to zero
//  std::fill(Corr.data(), Corr.data()+Corr.num_elements(), cmplx(.0,.0));
  std::fill(C4_mes.data(), C4_mes.data()+C4_mes.num_elements(), cmplx(.0,.0));

  // global variables from input file needed here
  const int Lt = global_data->get_Lt();

  // memory for intermediate matrices when building C4_3 (save multiplications)
  LapH::CrossOperator X(2);

  basic.init_operator('b', vdaggerv, peram);

  // computing the meson correlator which can be used to compute all small
  // trace combinations for 2pt and 4pt functions
  build_Corr();

  // computing the meson 4pt big cross trace
  // TODO: if condition that at least four random vectos are needed
  compute_meson_4pt_cross_trace(X);
//
  write_C4_3(config_i);
  build_and_write_2pt(config_i);
  build_and_write_C4_1(config_i);
  build_and_write_C4_2(config_i);

}

void grsource_w() {
register int i,j,k;
register site *s;
    FORMYSITES(i,s){
        for(k=0;k<4;k++)for(j=0;j<3;j++){
#ifdef SITERAND
            s->g_rand.d[k].c[j].real = gaussian_rand_no(&(s->site_prn));
            s->g_rand.d[k].c[j].imag = gaussian_rand_no(&(s->site_prn));
#else
            s->g_rand.d[k].c[j].real = gaussian_rand_no(&node_prn);
            s->g_rand.d[k].c[j].imag = gaussian_rand_no(&node_prn);
#endif
	    s->psi.d[k].c[j] = cmplx(0.0,0.0);
        }
    }

#ifdef LU
    /* use g_rand on odd sites as temporary storage:
        g_rand(odd) = Dslash^adjoint * g_rand(even) */
    dslash_w_site( F_OFFSET(g_rand), F_OFFSET(g_rand), MINUS, ODD);
    dslash_w_site( F_OFFSET(g_rand), F_OFFSET(chi)   , MINUS, EVEN);
    FOREVENSITES(i,s){
        scalar_mult_add_wvec( &(s->g_rand), &(s->chi), -kappa*kappa,
            &(s->chi) );
    }

//Processing when import frame
void frame_in(){
	int k;
	float mag;
	for(k=0;k<FFTLEN;k++){
		buffer[k]=cmplx(inframe[k],0);	
	}	 
	//perform FFT
	fft(FFTLEN,buffer);
	//rotate every 2.5s
	if(i== 312){
	M4=M3;
	M3=M2;
	M2=M1;	
	M1=M4;						      		
	i=0;
	}
	//Import to M1
	if(i==0){
		for(k=0;k<FFTLEN;k++){
			M1[k]=cabs(buffer[k]);			
		}	
	}
	else{
		for(k=0;k<FFTLEN;k++){
			mag=cabs(buffer[k]);
			M1[k]=(1-lpf_weight_speech)*mag+lpf_weight_speech*M1[k];
		}
	}
}

main()
{
	int n1,n2,j;
	float scale,sumz,sume;
 
	for (n1=1,n2=npfa(n1+1); n2<NMAX; n1=n2,n2=npfa(n1+1)) {
		for (j=0; j<n1*n2; j++)
			c[j] = z[j] = cmplx(franuni(),franuni());
		pfamcc(1,n1,n2,1,n1,z);
		pfamcc(1,n2,n1,n1,1,z);
		pfamcc(-1,n1,n2,1,n1,z);
		pfamcc(-1,n2,n1,n1,1,z);
		scale = 1.0/(float)(n1*n2);
		for (j=0; j<n1*n2; j++)
			z[j] = crmul(z[j],scale);

		for (j=0,sumz=sume=0.0; j<n1*n2; j++) {
			sumz += fcabs(z[j]);
			sume += fcabs(csub(z[j],c[j]));
		}
		printf("n1 = %d  n2 = %d  sume/sumz = %0.10f\n",
			n1,n2,sume/sumz);
		if (sume/sumz>1.0e-6) printf("!!! warning !!!\n");
	}
}

main()
{
	int n,nc,nr,i;
	float *rz=(float*)cz;
	float cpucc,cpurc;

	for (nc=npfa(NMIN),nr=npfar(nc);
		nc<NMAX && nr<NMAX;
		nc=npfa(nc+1),nr=npfar(nc)) {

		for (i=0; i<nc*nc; i++)
			cz[i] = cmplx(0.0,0.0);
		cpucc = cpusec();
		pfa2cc(1,1,nc,nc,cz);
		cpucc = cpusec()-cpucc;

		for (i=0; i<nr*nr; i++)
			rz[i] = 0.0;
		cpurc = cpusec();
		pfa2rc(1,1,nr,nr,rz,cz);
		cpurc = cpusec()-cpurc;

		printf("nc,nr,cc,rc,cc/rc = %d %d %f %f %f\n",
			nc,nr,cpucc,cpurc,cpucc/cpurc);
	}
}

// -----------------------------------------------------------------------------
// -----------------------------------------------------------------------------
void LapH::Correlators::build_C1(const Quarklines& quarklines,
                    const std::vector<CorrInfo>& corr_lookup,
                    const QuarklineLookup& quark_lookup,
                    const std::vector<RandomIndexCombinationsQ2>& ric_lookup) {

  std::cout << "\tcomputing C1:";
  clock_t time = clock();

  for(const auto& c_look : corr_lookup){
    const auto& ric = ric_lookup[quark_lookup.Q1[c_look.lookup[0]].
                                                     id_ric_lookup].rnd_vec_ids;
    std::vector<cmplx> correlator(Lt*ric.size(), cmplx(.0,.0));
    for(size_t t = 0; t < Lt; t++){
      for(const auto& id : ric){
        correlator[(&id-&ric[0])*Lt + t] +=
         quarklines.return_Q1(t, t/dilT, c_look.lookup[0], &id-&ric[0]).trace();
      }
    }
    // write data to file
    write_correlators(correlator, c_look);
  }
  time = clock() - time;
  std::cout << "\t\t\tSUCCESS - " << ((float) time) / CLOCKS_PER_SEC 
            << " seconds" << std::endl;
}

// -----------------------------------------------------------------------------
// -----------------------------------------------------------------------------
void LapH::Correlators::build_C2c(const std::vector<CorrInfo>& corr_lookup) {

  for(const auto& c_look : corr_lookup){
    std::vector<cmplx> correlator(Lt, cmplx(.0,.0));
    if(c_look.outfile.find("Check") == 0){
      for(int t1 = 0; t1 < Lt; t1++){
        for(const auto& corr : corrC[c_look.lookup[0]][t1][t1]){
          correlator[t1] += corr;
        }
      }
      // normalisation
      for(auto& corr : correlator)
        corr /= corrC[c_look.lookup[0]][0][0].size();
      // write data to file
      write_correlators(correlator, c_look);
    }
    else{
      for(int t1 = 0; t1 < Lt; t1++){
      for(int t2 = 0; t2 < Lt; t2++){
        int t = abs((t2 - t1 - (int)Lt) % (int)Lt);
        for(const auto& corr : corrC[c_look.lookup[0]][t1][t2]){
          correlator[t] += corr;
        }
      }}
      // normalisation
      for(auto& corr : correlator)
        corr /= Lt*corrC[c_look.lookup[0]][0][0].size();
      // write data to file
      write_correlators(correlator, c_look);
    }
  }
}

void rectangle::divide(rectangle & r1, rectangle & r2, double ratio){
    auto width = (z2 - z1).real();
    auto height = (z2 - z1).imag();

    // вот этот кусок стоит упростить
    if (width > height) {
        auto u1 = z1, u2 = z2 - (z2 - z1).real() / (1. + ratio);
        r1 = {u1, u2};

        u1 = z1 + (z2 - z1).real() * ratio / (1. + ratio), u2 = z2;
        r2 = {u1, u2};
    } else {
        auto u1 = z1, u2 = z2 - cmplx(0, (z2 - z1).imag() / (1. + ratio));
        r1 = {u1, u2};

        u1 = z1 + cmplx(0, (z2 - z1).imag() * ratio / (1. + ratio)), u2 = z2;
        r2 = {u1, u2};
    }
}

void wl_2l_2corr_offax(int tot_smear, int step) {

    register int i,dir1,dir2,dir3,r,t,r_off;
    int nth,nxh,nrmax,j,k,num_lp,rmax;
    register site *s;
    dble_su3_vec_src *dpt;
    su3_vector tvec1,tvec2;
    msg_tag *mtag[8],*gmtag;
    complex  cc,cc1,cc2;
    complex *wl2_mes, *wl2_mes_d;
    Real ftmpr, ftmpi, fns, xnsrc;
    field_offset meson1, meson2, meson_f, mes_t;
    int disp[4];    /* displacement vector for general gather */


    meson1 = F_OFFSET(xxx.c[0]);
    meson2 = F_OFFSET(xxx.c[1]);
    meson_f = F_OFFSET(ttt.c[0]);
    mes_t = F_OFFSET(ttt.c[1]);

    if( nx != ny || nx != nz) {
        if(this_node == 0)printf("wl_2lprop gives wrong results for nx!=ny!=nz");
        return;
    }

    xnsrc = num_src*(num_src-1);
    nth = nt/2;
    nxh = nx/2;
    nrmax = 4;
    wl2_mes = (complex *)malloc(nth*nrmax*sizeof(complex));
    wl2_mes_d = (complex *)malloc(nth*nrmax*sizeof(complex));

    for(t=0; t<nth; t++) for(r=0; r<nrmax; r++) {
            wl2_mes[r+nrmax*t] = cmplx(0.0,0.0);
            wl2_mes_d[r+nrmax*t] = cmplx(0.0,0.0);
        }


    FORALLSITES(i,s) {
        su3_vec_to_src( &(s->qprop[0]), &(s->dtmpvecs[0].n[0]), num_src);
        su3_vec_to_src( &(s->g_rand[0]), &(s->dtmpvecs[0].n[1]), num_src);
    }

void randomize(field_offset G, Real radius)
{
  register int a, b, i;
  site *s;

  FORALLSITES(i,s) {
    for(a=0; a<3; a++) for(b=0; b<3; b++)
      (*(su3_matrix *)F_PT(s,G)).e[a][b] 
             = cmplx(radius*((Real)drand48()-0.5),
                     radius*((Real)drand48()-0.5));
    reunit_su3((su3_matrix *)F_PT(s,G));
  }
}

static void 
get_Q_from_VUadj(su3_matrix *Q, su3_matrix *V, su3_matrix *U){

  complex minusI;
  complex tr;
  su3_matrix Omega;

  minusI = cmplx(0, -1);  /* -i */

  /* Omega = V U^adj */
  mult_su3_na( V, U, &Omega );

  /* Q = traceless hermitian part of [-i Omega] */
  c_scalar_mult_su3mat( &Omega, &minusI, &Omega );
  traceless_hermitian_su3( Q, &Omega );
}

//Processing when import frame
void frame_in(){
	int k,n;
	float mag;
	
	for(k=0;k<FFTLEN;k++){
		buffer[k]=cmplx(inframe[k],0);	
	}	 
	//perform FFT
	fft(FFTLEN,buffer);
	//rotate every 2.5s
	n=rotate_inteval*FSAMP*OVERSAMP/FFTLEN;
	if(i== n){
	M4=M3;
	M3=M2;
	M2=M1;	
	M1=M4;						      		
	i=0;
	}
	//Import to M1
	if(i==0){
		for(k=0;k<FFTLEN;k++){
			M1[k]=cabs(buffer[k]);			
		}	
	}
	else{
		if(lpf_switch==0){
			for(k=0;k<FFTLEN;k++){
				mag=cabs(buffer[k]);
				if(M1[k]>mag){
					M1[k]=mag;
				}
			}
		}
		else{
			for(k=0;k<FFTLEN;k++){
				if(lpf_power_switch==1){
				mag=cabs(cmul(buffer[k],buffer[k]));
				M1[k]=sqrt((1-lpf_weight)*mag+lpf_weight*(M1[k]*M1[k]));
				}
				else{
				mag=cabs(buffer[k]);
				M1[k]=(1-lpf_weight)*mag+lpf_weight*M1[k];
				}
			}
		}	 
	}
}

main()
{
	int n=NMIN,nfft,i,it;
	float cpucc;

	for (i=0; i<NMAX; i++)
		z[i] = cmplx(0.0,0.0);

	for (n=NMIN; n<=NMAX; n+=NSTEP) {
		nfft = npfa(n);
		cpucc = cpusec();
		for (it=0; it<NT; it++)
			pfacc(1,nfft,z);
		cpucc = cpusec()-cpucc;
		cpucc /= NT;
		printf("n = %d  nfft = %d  sec = %f\n",n,nfft,cpucc);
	}
}