#include<stdio.h>
#include<math.h>

/*****************************************************************************
*                             Program of Setup                               *
*****************************************************************************/
#define dx    (0.3)
#define dy    (0.3)
#define imf   (200)
#define jmf   (100)
#define PI    (3.141592)
#define dfx   (2.0*3.141592654/(dx*(imf-1)))
#define dfy   (2.0*3.141592654/(dy*(jmf-1)))
#define tg    (imf*jmf)

static double Re[imf][jmf],Xr[imf][jmf];
static double Im[imf][jmf],Xi[imf][jmf];
static double Tr[imf][jmf],Ti[imf][jmf];
static double fd[imf][jmf];

/*****************************************************************************
*                             Main program                                   *
*****************************************************************************/
int main(int argc,char *argv[])
{
  double pr,pi,tp1,tp2;
  int i,j,k;
  FILE *fp,*rp;

  pr = 2.0*3.141592/(double)(imf);
  pi = 2.0*3.141592/(double)(jmf);

  for(i=0;i<imf;i++)
    for(j=0;j<jmf;j++)
      {
	Re[i][j] = Im[i][j] = fd[i][j] = 0.0;
	Xr[i][j] = Xi[i][j] = Tr[i][j] = Ti[i][j] = 0.0;
      }

  if(argc != 3)
    {
      printf("Usage is ...\n");
      printf("dft.out <input file> <output file>\n");
      return -1;
    }

  /* reading of data file */
  if((rp = fopen(argv[1],"r")) == NULL) return -1;
  if((fp = fopen(argv[2],"w")) == NULL) return -1;

  printf("now loading ...\n");

  for(i=0;i<imf;i++)
    for(j=0;j<jmf;j++)
      fscanf(rp,"%lf %lf %lf\n",&tp1,&tp2,&fd[i][j]);

  printf("now DFT  process ...\n");

  /* 2 Dimention of DFT */
  for(j=0;j<jmf;j++)
    for(k=0;k<imf;k++)
      for(i=0;i<imf;i++)
	{
	  Tr[k][j] += fd[i][j]*cos(i*k*pr)/(double)imf;
	  Ti[k][j] -= fd[i][j]*sin(i*k*pr)/(double)imf;
	}
  
  for(k=0;k<imf;k++)
    for(i=0;i<jmf;i++)
      for(j=0;j<jmf;j++)
	{
	  Re[k][i] += Tr[k][j]*cos(i*j*pi) + Ti[k][j]*sin(i*j*pi);
	  Im[k][i] -= Tr[k][j]*sin(i*j*pi) - Ti[k][j]*cos(i*j*pi);
	}
  
   for(i=0;i<imf;i++)
     for(j=0;j<jmf;j++)
       {
	 Re[i][j] /= (double)jmf;
	 Im[i][j] /= (double)jmf;
       }

  printf("now IDFT process ...\n");
  
  for(i=0;i<imf;i++)
    for(j=0;j<jmf;j++)
      Tr[i][j] = Ti[i][j] = 0.0;
  
  /* 2 Dimention of IDFT */
  for(j=0;j<jmf;j++)
    for(k=0;k<imf;k++)
      for(i=0;i<imf;i++)
	{
	  Tr[k][j] += Re[i][j]*cos(i*k*pr) - Im[i][j]*sin(i*k*pr);
	  Ti[k][j] += Re[i][j]*sin(i*k*pr) + Im[i][j]*cos(i*k*pr);
	}
  
  for(k=0;k<imf;k++)
    for(i=0;i<jmf;i++)
      for(j=0;j<jmf;j++)
	{
	  Xr[k][i] += Tr[k][j]*cos(i*j*pi) - Ti[k][j]*sin(i*j*pi);
	  Xi[k][i] += Tr[k][j]*sin(i*j*pi) + Ti[k][j]*cos(i*j*pi);
	}

  for(i=0;i<imf;i++){
    for(j=0;j<jmf;j++)
      fprintf(fp,"%e %e %e %e %e %e %e %e\n",i*dx,j*dy,i*dfx,j*dfy,
	      fd[i][j],Re[i][j]*Re[i][j]+Im[i][j]*Im[i][j],Xr[i][j],Xi[i][j]);
    fprintf(fp,"\n");
  }
  
  printf("... complete\n");

  fclose(fp);
  fclose(rp);

  return 0;
}