// to show the Nyquist frequency 
// we sample the same signal twice, first above the critical
// frequency and later below
// f(t) = sin(2*pi*4*t)   f = 4Hz
//
// time interval: 2*pi sec
//  first sampling: 64 points, delta_t = 0.1 sec, f_c = 5 Hz
//  second sampling: 32 points, delta_t = 0.2 sec, f_c = 2.5 Hz
//      jtmn, Nov, 2003

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_math.h>

#define REAL(z,i) ((z)[2*(i)])
#define IMAG(z,i) ((z)[2*(i)+1])

int main (void)
{
  int i;
  double dataok[2*64];
  double dataprob[2*32];

  double t = 0.0;
  for (i = 0; i < 32; i++)
    {
      // delta_t = 0.2 -> Nyquist freq = 2.5 Hz
      REAL(dataprob,i) = sin(2*M_PI*4.0*t);
      IMAG(dataprob,i) = 0.0;
      t = t + 0.2;
    }

  t=0.0;
  for (i = 0; i < 64; i++)
    {
      // delta_t = 0.1 -> Nyquist freq = 5 Hz
      REAL(dataok,i) = sin(2*M_PI*4.0*t);
      IMAG(dataok,i) = 0.0;
      t = t + 0.1;
    }

  FILE *in_file1;
  in_file1 = fopen ("inputprob.dat", "w");
  for (i= 0; i < 32; i++)
    {
     fprintf (in_file1,"%d %e %e\n", i, REAL(dataprob,i), IMAG(dataprob,i));  
    }
  printf ("\n");

  FILE *in_file2;
  in_file2 = fopen ("inputok.dat", "w");
  for (i= 0; i < 64; i++)
    {
     fprintf (in_file2,"%d %e %e\n", i, REAL(dataok,i), IMAG(dataok,i));  
    }
  printf ("\n");

  gsl_fft_complex_radix2_forward (dataprob, 1, 32);
  gsl_fft_complex_radix2_forward (dataok, 1, 64);

  FILE *out_file1;
  out_file1 = fopen ("outputprob.dat", "w");
  for (i = 0; i < 32; i++)
    {
      fprintf (out_file1, "%d %e %e\n", i, 
	       REAL(dataprob,i),
	       IMAG(dataprob,i));
    }

  FILE *out_file2;
  out_file2 = fopen ("outputok.dat", "w");
  for (i = 0; i < 64; i++)
    {
      fprintf (out_file2, "%d %e %e\n", i, 
	       REAL(dataok,i),
	       IMAG(dataok,i));
    }  

  return 0; 
}