rfft.cpp

The Spectral Toolkit is a C++ spectral transform library written by Rodney James and Chuck Panaccion

	  z2.re=t0.re+a1.re+a2.im; 
	  z2.im=t0.im-a1.im-a2.re; 
	  x0[L-1]=2.0*t0.re+t1.re;
	  x1[L-1]=z1.re*w1.re+z1.im*w1.im;
	  x2[L-1]=z2.re*w2.re+z2.im*w2.im;
        }
    }
}

/// Internal radix-4 butterfly forward transform.
/// \param x input array
/// \param y output array
/// \param L accumulation product of previous factors
/// \param N remaining sequence length from node down

void rfft::node::analysis_radix4(real *x,real *y,int L,int N)
{
  complex t0,t1,t2,t3;
  complex w,w2,w3;
  complex z0,z1,z2,z3;
  real *x0,*x1,*x2,*x3,*y0,*y1,*y2;
    
  for(int j=0;j<N;j++)
    {
      x0=x+j*L;
      x1=x+j*L+N*L;
      x2=x+j*L+2*N*L;
      x3=x+j*L+3*N*L;
      y0=y+j*4*L;
      y1=y+j*4*L+2*L;
      y2=y+j*4*L+4*L;
      t0.re=x0[0];
      t1.re=x1[0];
      t2.re=x2[0];
      t3.re=x3[0];
      z1.re=t0.re+t2.re;
      z1.im=t1.re+t3.re;
      y0[ 0]= z1.re+z1.im;
      y1[-1]= t0.re-t2.re;
      y1[ 0]=-t1.re+t3.re;      
      y2[-1]= z1.re-z1.im;
      w.re=omega.re+1.0;                   
      w.im=omega.im;
      for(int k=1;k<(L+1)/2;k++)
        {
	  w2.re=w.re*w.re-w.im*w.im;
	  w2.im=2.0*w.re*w.im;
	  w3.re=w.re*w2.re-w.im*w2.im;
	  w3.im=w.re*w2.im+w.im*w2.re;
	  t0.re=x0[2*k-1];
	  t0.im=x0[2*k  ];
	  z1.re=x1[2*k-1];
	  z1.im=x1[2*k  ];
	  z2.re=x2[2*k-1];
	  z2.im=x2[2*k  ];
	  z3.re=x3[2*k-1];
	  z3.im=x3[2*k  ];
	  t1.re=z1.re*w.re-z1.im*w.im;
	  t1.im=z1.im*w.re+z1.re*w.im;
	  t2.re=z2.re*w2.re-z2.im*w2.im;
	  t2.im=z2.im*w2.re+z2.re*w2.im;
	  t3.re=z3.re*w3.re-z3.im*w3.im;
	  t3.im=z3.im*w3.re+z3.re*w3.im;                        
	  z0.re=t0.re+t2.re;
	  z0.im=t0.im+t2.im;
	  z1.re=t0.re-t2.re;
	  z1.im=t0.im-t2.im;
	  z2.re=t1.re+t3.re;
	  z2.im=t1.im+t3.im;
	  z3.re=t1.re-t3.re;
	  z3.im=t1.im-t3.im;
	  y0[ 2*k-1]= (z0.re+z2.re);
	  y0[ 2*k  ]= (z0.im+z2.im);
	  y1[ 2*k-1]= (z1.re+z3.im);
	  y1[ 2*k  ]= (z1.im-z3.re);
	  y1[-2*k-1]= (z1.re-z3.im);
	  y1[-2*k  ]=-(z1.im+z3.re);
	  y2[-2*k-1]= (z0.re-z2.re);
	  y2[-2*k  ]=-(z0.im-z2.im);
	  z1.re=omega.re*w.re-omega.im*w.im+w.re;
	  z1.im=omega.re*w.im+omega.im*w.re+w.im; 
	  w.re=z1.re;
	  w.im=z1.im;
        }
      if(L%2==0)
        {
	  w2.re=w.re*w.re-w.im*w.im;
	  w2.im=2.0*w.re*w.im;
	  w3.re=w.re*w2.re-w.im*w2.im;
	  w3.im=w.re*w2.im+w.im*w2.re;
	  t0.re=x0[L-1];
	  t0.im=0.0;
	  t1.re=x1[L-1]*w.re;
	  t1.im=x1[L-1]*w.im;
	  t2.re=x2[L-1]*w2.re;
	  t2.im=x2[L-1]*w2.im;
	  t3.re=x3[L-1]*w3.re;
	  t3.im=x3[L-1]*w3.im;
	  z0.re=t0.re+t2.re;
	  z0.im=t0.im+t2.im;
	  z1.re=t0.re-t2.re;
	  z1.im=t0.im-t2.im;
	  z2.re=t1.re+t3.re;
	  z2.im=t1.im+t3.im;
	  z3.re=t1.re-t3.re;
	  z3.im=t1.im-t3.im;
	  y0[L-1]=z0.re+z2.re;
	  y0[L  ]=z0.im+z2.im;
	  y1[L-1]=z1.re+z3.im;
	  y1[L  ]=z1.im-z3.re;
        }
    }
}

/// Internal radix-4 butterfly backward transform.
/// \param x input array
/// \param y output array
/// \param L accumulation product of previous factors
/// \param N remaining sequence length from node down

void rfft::node::synthesis_radix4(real *y,real *x,int L,int N)
{
  complex t0,t1,t2,t3;
  complex w,w2,w3;
  complex z0,z1,z2,z3;
  real *x0,*x1,*x2,*x3,*y0,*y1,*y2;
    
  for(int j=0;j<N;j++)
    {
      x0=x+j*L;
      x1=x+j*L+N*L;
      x2=x+j*L+2*N*L;
      x3=x+j*L+3*N*L;
      y0=y+j*4*L;
      y1=y+j*4*L+2*L;
      y2=y+j*4*L+4*L;
      t0.re=y0[ 0];
      t1.re=y1[-1];
      t1.im=y1[ 0];
      t2.re=y2[-1];
      z0.re=t0.re+t2.re;
      z0.im=t0.re-t2.re;
      z1.re=2.0*t1.re;
      z1.im=2.0*t1.im;
      x0[0]=z0.re+z1.re;
      x1[0]=z0.im-z1.im; 
      x2[0]=z0.re-z1.re;
      x3[0]=z0.im+z1.im;
      w.re=omega.re+1.0;                   
      w.im=omega.im;
      for(int k=1;k<(L+1)/2;k++)
        {      
	  w2.re=w.re*w.re-w.im*w.im;
	  w2.im=2.0*w.re*w.im;
	  w3.re=w.re*w2.re-w.im*w2.im;
	  w3.im=w.re*w2.im+w.im*w2.re;
	  t0.re= y0[ 2*k-1];
	  t0.im= y0[ 2*k  ];
	  t1.re= y1[ 2*k-1];
	  t1.im= y1[ 2*k  ];
	  t2.re= y1[-2*k-1];
	  t2.im=-y1[-2*k  ];
	  t3.re= y2[-2*k-1];
	  t3.im=-y2[-2*k  ];
	  z0.re=t0.re+t3.re;
	  z0.im=t0.im+t3.im;
	  z1.re=t0.re-t3.re;
	  z1.im=t0.im-t3.im;
	  z2.re=t1.re+t2.re;
	  z2.im=t1.im+t2.im;
	  z3.re=t1.re-t2.re;
	  z3.im=t1.im-t2.im;
	  t0.re=z0.re+z2.re;
	  t0.im=z0.im+z2.im;
	  t1.re=z1.re-z3.im;
	  t1.im=z1.im+z3.re;
	  t2.re=z0.re-z2.re;
	  t2.im=z0.im-z2.im;
	  t3.re=z1.re+z3.im;
	  t3.im=z1.im-z3.re;
	  x0[2*k-1]=t0.re;
	  x0[2*k  ]=t0.im;
	  x1[2*k-1]=t1.re*w.re+t1.im*w.im;
	  x1[2*k  ]=t1.im*w.re-t1.re*w.im;
	  x2[2*k-1]=t2.re*w2.re+t2.im*w2.im;
	  x2[2*k  ]=t2.im*w2.re-t2.re*w2.im;
	  x3[2*k-1]=t3.re*w3.re+t3.im*w3.im;
	  x3[2*k  ]=t3.im*w3.re-t3.re*w3.im;
	  z0.re=omega.re*w.re-omega.im*w.im+w.re;
	  z0.im=omega.re*w.im+omega.im*w.re+w.im; 
	  w.re=z0.re;
	  w.im=z0.im;
        }
      if(L%2==0)
        {
	  w2.re=w.re*w.re-w.im*w.im;
	  w2.im=2.0*w.re*w.im;
	  w3.re=w.re*w2.re-w.im*w2.im;
	  w3.im=w.re*w2.im+w.im*w2.re;
	  t0.re=y0[L-1];
	  t0.im=y0[L  ];
	  t1.re=y1[L-1];
	  t1.im=y1[L  ];
	  t2.re= t1.re;
	  t2.im=-t1.im;
	  t3.re= t0.re;
	  t3.im=-t0.im;         
	  z0.re=t0.re+t2.re;
	  z0.im=t0.im+t2.im;
	  z1.re=t0.re-t2.re;
	  z1.im=t0.im-t2.im;
	  z2.re=t1.re+t3.re;
	  z2.im=t1.im+t3.im;
	  z3.re=t1.re-t3.re;
	  z3.im=t1.im-t3.im;
	  t0.re=z0.re+z2.re;
	  t1.re=z1.re-z3.im;
	  t1.im=z1.im+z3.re;
	  t2.re=z0.re-z2.re;
	  t2.im=z0.im-z2.im;
	  t3.re=z1.re+z3.im;
	  t3.im=z1.im-z3.re;
	  x0[L-1]=t0.re;
	  x1[L-1]=t1.re*w.re+t1.im*w.im;
	  x2[L-1]=t2.re*w2.re+t2.im*w2.im;
	  x3[L-1]=t3.re*w3.re+t3.im*w3.im;
        }
    }
}

/// Internal radix-5 butterfly forward transform.
/// \param x input array
/// \param y output array
/// \param L accumulation product of previous factors
/// \param N remaining sequence length from node down

void rfft::node::analysis_radix5(real *x,real *y,int L,int N)
{
  complex t0,t1,t2,t3,t4;
  complex w1,w2,w3,w4;
  complex z0,z1,z2,z3,z4;
  real *x0,*x1,*x2,*x3,*x4,*y0,*y1,*y2;
  const real a1=0.3090169943749474241022934171830L;
  const real b1=-0.9510565162951535721164393333793L;
  const real a2=-0.8090169943749474241022934171825L;
  const real b2=-0.5877852522924731291687059546394L;
    
  for(int j=0;j<N;j++)
    {
      x0=x+j*L;
      x1=x+j*L+N*L;
      x2=x+j*L+2*N*L;
      x3=x+j*L+3*N*L;
      x4=x+j*L+4*N*L;
      y0=y+j*5*L;
      y1=y+j*5*L+2*L;
      y2=y+j*5*L+4*L;
        
      t0.re=x0[0];
      t1.re=x1[0];
      t2.re=x2[0];
      t3.re=x3[0];
      t4.re=x4[0];
      z1.re=t1.re+t4.re;
      z2.re=t1.re-t4.re;
      z3.re=t2.re+t3.re;
      z4.re=t2.re-t3.re;
        
      y0[ 0]=t0.re+z1.re+z3.re;      
      y1[-1]=t0.re+a1*z1.re+a2*z3.re;
      y1[ 0]=      b1*z2.re+b2*z4.re;  
      y2[-1]=t0.re+a2*z1.re+a1*z3.re;
      y2[ 0]=      b2*z2.re-b1*z4.re;
      w1.re=omega.re+1.0;                   
      w1.im=omega.im;
      for(int k=1;k<(L+1)/2;k++)
        {
	  w2.re = w1.re*w1.re-w1.im*w1.im;
	  w2.im = 2.0*w1.re*w1.im;
	  w3.re = w2.re*w1.re-w2.im*w1.im;
	  w3.im = w2.re*w1.im+w2.im*w1.re; 
	  w4.re = w3.re*w1.re-w3.im*w1.im;
	  w4.im = w3.re*w1.im+w3.im*w1.re;
	  z0.re=x0[2*k-1];
	  z0.im=x0[2*k  ];
	  z1.re=x1[2*k-1];
	  z1.im=x1[2*k  ];
	  z2.re=x2[2*k-1];
	  z2.im=x2[2*k  ];
	  z3.re=x3[2*k-1];
	  z3.im=x3[2*k  ];
	  z4.re=x4[2*k-1];
	  z4.im=x4[2*k  ];
	  t1.re=z1.re*w1.re-z1.im*w1.im;
	  t1.im=z1.im*w1.re+z1.re*w1.im;
	  t2.re=z2.re*w2.re-z2.im*w2.im;
	  t2.im=z2.im*w2.re+z2.re*w2.im;
	  t3.re=z3.re*w3.re-z3.im*w3.im;
	  t3.im=z3.im*w3.re+z3.re*w3.im;
	  t4.re=z4.re*w4.re-z4.im*w4.im;
	  t4.im=z4.im*w4.re+z4.re*w4.im;
	  z1.re=t1.re+t4.re;
	  z1.im=t1.im+t4.im;
	  z2.re=t1.re-t4.re;
	  z2.im=t1.im-t4.im;
	  z3.re=t2.re+t3.re;
	  z3.im=t2.im+t3.im;
	  z4.re=t2.re-t3.re;
	  z4.im=t2.im-t3.im;
	  t0.re=z0.re+z1.re+z3.re;
	  t0.im=z0.im+z1.im+z3.im;
	  t1.re=z0.re+a1*z1.re+a2*z3.re;
	  t1.im=z0.im+a1*z1.im+a2*z3.im;
	  t2.re=z0.re+a1*z3.re+a2*z1.re;
	  t2.im=z0.im+a1*z3.im+a2*z1.im;
	  t3.re=b1*z2.re+b2*z4.re;  
	  t3.im=b1*z2.im+b2*z4.im;
	  t4.re=b1*z4.re-b2*z2.re;
	  t4.im=b1*z4.im-b2*z2.im;
	  y0[ 2*k-1]= t0.re;
	  y0[ 2*k  ]= t0.im;
	  y1[ 2*k-1]= t1.re-t3.im;
	  y1[ 2*k  ]= t1.im+t3.re;
	  y2[ 2*k-1]= t2.re+t4.im;
	  y2[ 2*k  ]= t2.im-t4.re;
	  y1[-2*k-1]= t1.re+t3.im;
	  y1[-2*k  ]=-t1.im+t3.re;
	  y2[-2*k-1]= t2.re-t4.im;
	  y2[-2*k  ]=-t2.im-t4.re;
	  z1.re=omega.re*w1.re-omega.im*w1.im+w1.re;
	  z1.im=omega.re*w1.im+omega.im*w1.re+w1.im; 
	  w1.re=z1.re;
	  w1.im=z1.im;
        }
      if(L%2==0)
        {
	  w2.re = w1.re*w1.re-w1.im*w1.im;
	  w2.im = 2.0*w1.re*w1.im;
	  w3.re = w2.re*w1.re-w2.im*w1.im;
	  w3.im = w2.re*w1.im+w2.im*w1.re; 
	  w4.re = w3.re*w1.re-w3.im*w1.im;
	  w4.im = w3.re*w1.im+w3.im*w1.re;
	  t0.re=x0[L-1];
	  t1.re=x1[L-1]*w1.re;
	  t1.im=x1[L-1]*w1.im;
	  t2.re=x2[L-1]*w2.re;
	  t2.im=x2[L-1]*w2.im;
	  t3.re=x3[L-1]*w3.re;
	  t3.im=x3[L-1]*w3.im;
	  t4.re=x4[L-1]*w4.re;
	  t4.im=x4[L-1]*w4.im;
	  z1.re=t1.re+t4.re;
	  z1.im=t1.im+t4.im;
	  z2.re=t1.re-t4.re;
	  z2.im=t1.im-t4.im;
	  z3.re=t2.re+t3.re;
	  z3.im=t2.im+t3.im;
	  z4.re=t2.re-t3.re;
	  z4.im=t2.im-t3.im;
	  y0[L-1]=t0.re+z1.re+z3.re;
	  y0[L  ]=z1.im+z3.im;
	  y1[L-1]=(t0.re+a1*z1.re-b1*z2.im+a2*z3.re-b2*z4.im);
	  y1[L  ]=(      b1*z2.re+a1*z1.im+b2*z4.re+a2*z3.im);
	  y2[L-1]=(t0.re+a2*z1.re-b2*z2.im+a1*z3.re+b1*z4.im);
        }
    }
}

/// Internal radix-5 butterfly backward transform.
/// \param x input array
/// \param y output array
/// \param L accumulation product of previous factors
/// \param N remaining sequence length from node down

void rfft::node::synthesis_radix5(real *y,real *x,int L,int N)
{
  complex t0,t1,t2,t3,t4;
  complex w1,w2,w3,w4;
  complex z0,z1,z2,z3,z4;
  real *x0,*x1,*x2,*x3,*x4,*y0,*y1,*y2;
  const real a1=0.3090169943749474241022934171830L;
  const real b1=-0.9510565162951535721164393333793L;
  const real a2=-0.8090169943749474241022934171825L;
  const real b2=-0.5877852522924731291687059546394L;
    
  for(int j=0;j<N;j++)
    {
      x0=x+j*L;
      x1=x+j*L+N*L;
      x2=x+j*L+2*N*L;
      x3=x+j*L+3*N*L;
      x4=x+j*L+4*N*L;
      y0=y+j*5*L;
      y1=y+j*5*L+2*L;
      y2=y+j*5*L+4*L;
        
      t0.re=y0[ 0];
      t1.re=2.0*y1[-1];
      t1.im=2.0*y1[ 0];
      t2.re=2.0*y2[-1];
      t2.im=2.0*y2[ 0];
      z1.re=(a1*t1.re+a2*t2.re);
      z1.im=(b1*t1.im+b2*t2.im);
      z2.re=(a2*t1.re+a1*t2.re);
      z2.im=(b2*t1.im-b1*t2.im);
      x0[0]=t0.re+(t1.re+t2.re);
      x1[0]=t0.re+(z1.re+z1.im);
      x2[0]=t0.re+(z2.re+z2.im);
      x3[0]=t0.re+(z2.re-z2.im);
      x4[0]=t0.re+(z1.re-z1.im);
      w1.re=omega.re+1.0;                   
      w1.im=omega.im;
      for(int k=1;k<(L+1)/2;k++)
        {
	  w2.re = w1.re*w1.re-w1.im*w1.im;
	  w2.im = 2.0*w1.re*w1.im;
	  w3.re = w2.re*w1.re-w2.im*w1.im;
	  w3.im = w2.re*w1.im+w2.im*w1.re; 
	  w4.re = w3.re*w1.re-w3.im*w1.im;
	  w4.im = w3.re*w1.im+w3.im*w1.re;
	  t0.re= y0[ 2*k-1];
	  t0.im= y0[ 2*k  ];
	  t1.re= y1[ 2*k-1];
	  t1.im= y1[ 2*k  ];
	  t2.re= y2[ 2*k-1];
	  t2.im= y2[ 2*k  ];
	  t3.re= y1[-2*k-1];
	  t3.im=-y1[-2*k  ];
	  t4.re= y2[-2*k-1];
	  t4.im=-y2[-2*k  ];
	  z0.re=t0.re;
	  z0.im=t0.im;
	  z1.re=t1.re+t3.re;
	  z1.im=t1.im+t3.im;
	  z2.re=t1.re-t3.re;
	  z2.im=t1.im-t3.im;
	  z3.re=t2.re+t4.re;
	  z3.im=t2.im+t4.im;
	  z4.re=t2.re-t4.re;
	  z4.im=t2.im-t4.im;
	  t0.re=z0.re+z1.re+z3.re;
	  t0.im=z0.im+z1.im+z3.im;
	  t1.re=z0.re+a1*z1.re+a2*z3.re; 
	  t1.im=z0.im+a1*z1.im+a2*z3.im;  
	  t2.re=z0.re+a1*z3.re+a2*z1.re;  
	  t2.im=z0.im+a1*z3.im+a2*z1.im;  
	  t3.re=b1*z2.re+b2*z4.re;
	  t3.im=b1*z2.im+b2*z4.im;
	  t4.re=b1*z4.re-b2*z2.re;
	  t4.im=b1*z4.im-b2*z2.im;
	  z0.re=t0.re;
	  z0.im=t0.im;
	  z1.re=t1.re+t3.im;
	  z1.im=t1.im-t3.re;
	  z2.re=t2.re-t4.im;
	  z2.im=t2.im+t4.re;
	  z3.re=t2.re+t4.im;
	  z3.im=t2.im-t4.re;
	  z4.re=t1.re-t3.im;
	  z4.im=t1.im+t3.re;
	  x0[2*k-1]=z0.re;
	  x0[2*k  ]=z0.im;
	  x1[2*k-1]=z1.re*w1.re+z1.im*w1.im;
	  x1[2*k  ]=z1.im*w1.re-z1.re*w1.im;
	  x2[2*k-1]=z2.re*w2.re+z2.im*w2.im;
	  x2[2*k  ]=z2.im*w2.re-z2.re*w2.im;      
	  x3[2*k-1]=z3.re*w3.re+z3.im*w3.im;
	  x3[2*k  ]=z3.im*w3.re-z3.re*w3.im;      
	  x4[2*k-1]=z4.re*w4.re+z4.im*w4.im;
	  x4[2*k  ]=z4.im*w4.re-z4.re*w4.im;
	  z1.re=omega.re*w1.re-omega.im*w1.im+w1.re;
	  z1.im=omega.re*w1.im+omega.im*w1.re+w1.im; 
	  w1.re=z1.re;
	  w1.im=z1.im;
        }
      if(L%2==0)
        {
	  w2.re = w1.re*w1.re-w1.im*w1.im;
	  w2.im = 2.0*w1.re*w1.im;
	  w3.re = w2.re*w1.re-w2.im*w1.im;
	  w3.im = w2.re*w1.im+w2.im*w1.re; 
	  w4.re = w3.re*w1.re-w3.im*w1.im;
	  w4.im = w3.re*w1.im+w3.im*w1.re;
	  t0.re=y0[L-1]; 
	  t0.im=y0[L  ]; 
	  t1.re=y1[L-1]; 
	  t1.im=y1[L  ]; 
	  t2.re=y2[L-1]; 
	  t2.im=0.0;
	  t3.re= t1.re;
	  t3.im=-t1.im;
	  t4.re= t0.re;
	  t4.im=-t0.im;
	  z0.re=t0.re;
	  z0.im=t0.im;
	  z1.re=t1.re+t4.re;
	  z1.im=t1.im+t4.im;
	  z2.re=t2.re+t3.re;
	  z2.im=t2.im+t3.im;
	  z3.re=t1.re-t4.re;
	  z3.im=t1.im-t4.im;
	  z4.re=t2.re-t3.re;
	  z4.im=t2.im-t3.im;
	  t0.re=z0.re+z1.re+z2.re;
	  t1.re=z0.re+a1*z1.re+a2*z2.re;
	  t1.im=z0.im+a1*z1.im+a2*z2.im;
	  t2.re=z0.re+a1*z2.re+a2*z1.re;
	  t2.im=z0.im+a1*z2.im+a2*z1.im;
	  t3.re=b1*z3.re+b2*z4.re;
	  t3.im=b1*z3.im+b2*z4.im;
	  t4.re=b1*z4.re-b2*z3.re;
	  t4.im=b1*z4.im-b2*z3.im;
	  z0.re=t0.re;
	  z1.re=t1.re+t3.im;
	  z1.im=t1.im-t3.re;
	  z2.re=t2.re-t4.im;
	  z2.im=t2.im+t4.re;
	  z3.re=t2.re+t4.im;
	  z3.im=t2.im-t4.re;
	  z4.re=t1.re-t3.im;
	  z4.im=t1.im+t3.re;
	  x0[L-1]=z0.re;    
	  x1[L-1]=z1.re*w1.re+z1.im*w1.im;
	  x2[L-1]=z2.re*w2.re+z2.im*w2.im;
	  x3[L-1]=z3.re*w3.re+z3.im*w3.im;
	  x4[L-1]=z4.re*w4.re+z4.im*w4.im;
        }
    }
}

/// Internal radix-6 butterfly forward transform.
/// \param x input array
/// \param y output array
/// \param L accumulation product of previous factors
/// \param N remaining sequence length from node down

void rfft::node::analysis_radix6(real *x,real *y,int L,int N)
{
  complex z,w,w2,w3,w4,w5;
  complex s0,s1,s2,s3,s4,s5;
  complex t0,t1,t2,t3,t4,t5;
  const real A= 0.5L;
  const real B=-0.5L*root3;
  real *x0,*x1,*x2,*x3,*x4,*x5;
  real *y0,*y1,*y2,*y3;
    
  for(int j=0;j<N;j++)
    {
      x0=x+j*L;
      x1=x+j*L+N*L;
      x2=x+j*L+2*N*L;
      x3=x+j*L+3*N*L;
      x4=x+j*L+4*N*L;
      x5=x+j*L+5*N*L;
      y0=y+j*6*L;
      y1=y+j*6*L+2*L;
      y2=y+j*6*L+4*L;
      y3=y+j*6*L+6*L;
      t0.re=x0[0];
      t1.re=x1[0];
      t2.re=x2[0];
      t3.re=x3[0];
      t4.re=x4[0];