/*  ------------------ apc.icc -------------------

   Example illustrating C programming for performing 
   operations comparable to those performed in apcxx.icc

                                   R. Ansari  08/2001
 
   ---- Computation steps :
    > Allocate an array representing a matrix (NL x NC )  (mtx)
    > fill it with a gaussian distributed random values
    > make a copy of the matrix   (mtxs)
    > allocate and fill 1D filter in Fourier space 
    > Loop over matrix rows k
     >> Extract row k (fline)
     >> compute 1D Fourier transform FFTForward (fline) 
     >> apply filter in Fourier space
     >> compute backward 1D FFT 
     >> Replace matrix row with the filtered values
   
    > free the allocated memory
   
   this example code can be 
     - included in a main program
     - executed using runcxx 
       csh> runcxx -tmpdir /tmp -f apc.icc
     - executed within spiapp 
       Cmd> c++execfrf apc.icc
*/

// Select computation on float or double (r_4 r_8)
#define FTYP r_4

// Number of matrix lines and colums
int NL = 1024;
int NC = 4096;
int i,k;
cout << " apc : NL= " << NL << " NC= " << NC << endl;
PrtTim("apc_Start");

// Creation of the initial matrix
FTYP * mtx = new FTYP [NL*NC];
// Filling matrix with gaussian random values
for(i=0; i<NL*NC; i++) 
  mtx[i] = 15.+NorRand()*3.;

// Allocation of the copy matrix (mtxs)
FTYP * mtxs = new FTYP [NL*NC];
// copying mtxs to mtx 
for(i=0; i<NL*NC; i++) 
  mtxs[i] = mtx[i];

// Creation and initialization of the Fourier filter filt(nu)
FTYP * filt = new FTYP [NC];
for(i=0; i<NC; i++) { 
  filt[i] =  1.-i/(double)NC;
  filt[i] *= filt[i]; 
}


// Creation of the FFTPackServer
// We use this class to avoid to do the initialization steps needed to use fftpack.
FFTPackServer ffts;
ffts.setNormalize(true);

PrtTim("apc_AfterInit");

// Allocation of the matrix row memory  
FTYP * fline = new FTYP [NC];

//  Loop over matrix rows
for(k=0; k<NL; k++) {
  // Matrix row extraction 
  for(i=0; i<NC; i++) fline[i] = mtx[k*NC+i];
  // Compute 1D forward FFT
  ffts.fftf(NC, fline);
  // Applying filter in Fourier space f(nu) =  f(nu)*filter(nu)
  for(i=0; i<NC; i++) fline[i] *= filt[i];
  // backward FFT
  ffts.fftb(NC, fline);
  // replace matrix row with filtered values
  for(i=0; i<NC; i++) mtx[k*NC+i] = fline[i];
}

PrtTim("apc_AfterFFTLoop");


delete[] mtx;
delete[] mtxs;
delete[] fline;
delete[] filt;