#include <iostream.h>
#include <math.h>
#include <iomanip.h>
#include <stdlib.h>
#include <fstream.h>
#include <string>
#include <process.h>
#include<stdio.h>

void prntft(double data[], double nn2)
{
	int n,m,mm;
    cout<<"n  Real H(n)  Imag H(n)  Real H(N-n)  Imag H(N-n)";
	cout<<endl;
    cout<<setw(1)<<"0";
	cout<<setw(14)<<data[1];
	cout<<setw(14)<<data[2];
	cout<<setw(14)<<data[1];
	cout<<setw(14)<<data[2]<<endl;
    for (n = 3; n<=(nn2 / 2) + 1; n=n+2)
	{
        m = (n - 1) / 2;
        mm = nn2 + 2 - n;
		coutsetiosflags(ios::fixed)<<;
        cout<<setprecision(0)<<setw(1)<<m;
		cout<<setprecision(6)<<setw(14)<<data[n];
		cout<<setprecision(6)<<setw(14)<<data[n+1];
		cout<<setprecision(6)<<setw(14)<<data[mm];
		cout<<setprecision(6)<<setw(14)<<data[mm+1]<<endl;
    }
}

void main()
{
    //program d12r1
    //driver for routine four1
	int nn,nn2,i,isign;
    double data[65], dcmp[65],aaa,j;
    nn = 32;
    nn2 = 2 * nn;   
    cout<<setw(1)<<"h(t)=real-valued even-function"<<endl;
    cout<<setw(1)<<"H(n)=H(N-n) and real?"<<endl;
    for (i = 1; i<=2 * nn - 1; i+=2)
	{
        data[i]=1.0/(((i-nn-1.0)/nn)*((i-nn-1.0)/nn)+1.0);
        data[i + 1] = 0.0;
    }
    isign = 1;
    four1(data, nn, isign);
    prntft(data, nn2);
    cout<<setw(1)<<"h(t)=imagnary-valued even-function"<<endl;
    cout<<setw(1)<<"H(n)=H(N-n) and imaginary?"<<endl;
    for (i = 1; i<= 2 * nn - 1; i+=2)
	{
        data[i+1]=1.0/(((i-nn-1.0)/nn)*((i-nn-1.0)/nn)+1.0);
        data[i] = 0.0;
    }
    isign = 1;
    four1(data, nn, isign);
    prntft(data, nn2);
	cout<<setw(1)<<"h(t)=real-valued odd-function"<<endl;
    cout<<setw(1)<<"H(n)=-H(N-n) and imaginary?"<<endl;
    for (i = 1; i<=2 * nn - 1; i+=2)
	{
        data[i]=(i-nn-1.0)/nn/(((i-nn-1.0)/nn)*((i-nn-1.0)/nn)+1.0);
        data[i + 1] = 0.0;
    }
    data[1] = 0.0;
    isign = 1;
    four1(data, nn, isign);
    prntft(data, nn2);
	cout<<setw(1)<<"h(t)=imagnary-valued odd-function"<<endl;
    cout<<setw(1)<<"H(n)=-H(N-n) and real?"<<endl;
    for( i = 1; i<= 2 * nn - 1; i+=2)
	{
        aaa = (i - nn - 1.0);
        data[i+1]=aaa/nn/(((i-nn-1.0)/nn)*((i-nn-1.0)/nn)+1.0);
        data[i] = 0.0;
    }
    data[2] = 0.0;
    isign = 1;
    four1(data, nn, isign);
    prntft(data, nn2);
    //transrorm, inverse-transform test
    for (i = 1; i<=2 * nn - 1; i+=2)
	{
        data[i]=1.0/((0.5*(i-nn-1)/nn)*(0.5*(i-nn-1)/nn)+1.0);
        dcmp[i] = data[i];
        data[i+1]=(0.25 * (i - nn - 1) / nn);
        data[i+1]=data[i+1]*exp(-(0.5*(i-nn-1)/nn)*(0.5*(i-nn-1)/nn));
        dcmp[i+1]=data[i + 1];
    }
    isign = 1;
    four1(data, nn, isign);
    isign = -1;
    four1(data, nn, isign);
	cout<<endl;
	cout<<" Double Fourier transform:      Original data:"<<endl;
    cout<<"k   Real h(k)   Imag h(k)   Real h(k)  Imag h(k)"<<endl;
    for (i = 1; i<=nn; i+=2)
	{
        j = (i + 1) / 2;
        cout<<setprecision(0)<<setw(1)<<j;
        cout<<setprecision(6)<<setw(14)<<dcmp[i];
        cout<<setprecision(6)<<setw(14)<<dcmp[i + 1];
        cout<<setprecision(6)<<setw(14)<<data[i]/nn;
        cout<<setprecision(6)<<setw(14)<<data[i + 1] / nn<<endl;
    }
}