#include <math.h>
#include "nr.h"
//--------------------------------------------------------------------
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
void fourn(double data[], unsigned long nn[], int ndim, int isign)
{int idim;
	unsigned long i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2;
	unsigned long ibit,k1,k2,n,nprev,nrem,ntot;
	double tempi,tempr;
	double theta,wi,wpi,wpr,wr,wtemp;

	for (ntot=1,idim=1;idim<=ndim;idim++)
		ntot *= nn[idim];
	nprev=1;
	for (idim=ndim;idim>=1;idim--) {
		n=nn[idim];
		nrem=ntot/(n*nprev);
		ip1=nprev << 1;
		ip2=ip1*n;
		ip3=ip2*nrem;
		i2rev=1;
		for (i2=1;i2<=ip2;i2+=ip1) {
			if (i2 < i2rev) {
				for (i1=i2;i1<=i2+ip1-2;i1+=2) {
					for (i3=i1;i3<=ip3;i3+=ip2) {
						i3rev=i2rev+i3-i2;
						SWAP(data[i3],data[i3rev]);
						SWAP(data[i3+1],data[i3rev+1]);
					}
				}
			}
			ibit=ip2 >> 1;
			while (ibit >= ip1 && i2rev > ibit) {
				i2rev -= ibit;
				ibit >>= 1;
			}
			i2rev += ibit;
		}
		ifp1=ip1;
		while (ifp1 < ip2) {
			ifp2=ifp1 << 1;
			theta=isign*6.28318530717959/(ifp2/ip1);
			wtemp=sin(0.5*theta);
			wpr = -2.0*wtemp*wtemp;
			wpi=sin(theta);
			wr=1.0;
			wi=0.0;
			for (i3=1;i3<=ifp1;i3+=ip1) {
				for (i1=i3;i1<=i3+ip1-2;i1+=2) {
					for (i2=i1;i2<=ip3;i2+=ifp2) {
						k1=i2;
						k2=k1+ifp1;
						tempr=(double)wr*data[k2]-(double)wi*data[k2+1];
						tempi=(double)wr*data[k2+1]+(double)wi*data[k2];
						data[k2]=data[k1]-tempr;
						data[k2+1]=data[k1+1]-tempi;
						data[k1] += tempr;
						data[k1+1] += tempi;
					}
				}
				wr=(wtemp=wr)*wpr-wi*wpi+wr;
				wi=wi*wpr+wtemp*wpi+wi;
			}
			ifp1=ifp2;
		}
		nprev *= n;
	}
}
#undef SWAP
/* (C) Copr. 1986-92 Numerical Recipes Software +-+. */
//--------------------------------------------------------------------


//--------------------------------------------------------------------
#define SWAP(a, b) tempr = (a); (a) = (b); (b) = tempr
void  four1(double data[], int nn, int isign)
{
	double tempr;
	int m;

	int n = nn << 1;
	int j = 1;
	for (int i = 1; i < n; i += 2) {
      		if (j > i) {
	  		SWAP(data[j],data[i]);
	  		SWAP(data[j+1],data[i+1]);
		};
      		m = n >> 1;
      		while (m >= 2 && j > m) {
	  		j  -= m; m >>= 1;
		};
      		j += m;
    	};
	int mmax = 2;
  	while (n > mmax) {
		double theta = (2*M_PI)/(isign*mmax);
		double wtemp = sin(0.5*theta);
		double wpr = -2.0*wtemp*wtemp;
		double wpi = sin(theta);
		double wr = 1.0;
		double wi = 0.0;
      		int istep = 2*mmax;
      		for (m = 1; m < mmax; m += 2) {
	  		for (int i = m; i <= n; i += istep) {
	      			j = i + mmax;
	      			tempr = wr*data[j ] - wi*data[j+1];
	      			double tempi = wr*data[j+1] + wi*data[j];
	      			data[j] = data[i] - tempr;
	      			data[j+1] = data[i+1] - tempi;
	      			data[i]	+= tempr;
	      			data[i+1] += tempi;
	    		};
	  		wr = (wtemp = wr)*wpr - wi*wpi + wr;
	  		wi = wi*wpr + wtemp*wpi + wi;
		};
		mmax = istep;
    	};
}
#undef SWAP
//--------------------------------------------------------------------
void realft(double data[], int n, int isign)
{
	double c1 = 0.5, c2, h1r, h1i, h2r, h2i;

	double theta = M_PI/(double)n;
	if (isign == 1) {
		c2 = -0.5; four1(data, n, 1);
	} else {
      		c2 =  0.5; theta = -theta;
	};
	double wtemp = sin(0.5*theta);
	double wpr = -2.0*wtemp*wtemp;
	double wpi = sin(theta);
	double wr = 1.0 + wpr;
	double wi = wpi;
	int n2p3 = 2*n + 3;
  	for (int i = 2; i <= n/2; i++) {
		int i1, i2, i3, i4;
		i4 =   1 + (i3 = n2p3 - (i2 = 1 + (i1 = i + i - 1)));
		h1r =  c1*(data[i1] + data[i3]);
		h1i =  c1*(data[i2] - data[i4]);
		h2r = -c2*(data[i2] + data[i4]);
		h2i =  c2*(data[i1] - data[i3]);
		data[i1] = h1r + wr*h2r - wi*h2i;
		data[i2] = h1i + wr*h2i + wi*h2r;
		data[i3] = h1r - wr*h2r + wi*h2i;
		data[i4] = -h1i + wr*h2i + wi*h2r;
		wr = (wtemp = wr)*wpr - wi*wpi + wr;
		wi = wi*wpr + wtemp*wpi + wi;
	}
	if (isign == 1) {
		data[1] = (h1r = data[1]) + data[2];
		data[2] = -data[2];
	} else {
		data[1] = c1*((h1r = data[1]) + data[2]);
		data[2] = c1*(h1r - data[2]);
		four1(data, n, -1);
	};
}
//--------------------------------------------------------------------
void cosft(double y[], int n, int isign)
{
	double enf0, even, odd, sume, sumo, y1, y2;
	double wi = 0.0;
	double wr = 1.0;
	double theta = M_PI/(double) n;
	double wtemp = sin(0.5*theta);
	double wpr = -2.0*wtemp*wtemp;
	double wpi = sin(theta);
	double sum = y[1];
	int m = n >> 1;
	int n2 = n + 2;
	int j;
	for (j = 2; j <= m; j++) {
		wr = (wtemp = wr)*wpr - wi*wpi + wr;
  		wi = wi*wpr + wtemp*wpi + wi;
		y1 = 0.5*(y[j] + y[n2-j]);
		y2 = (y[j] - y[n2-j]);
		y[j] = y1 - wi*y2;
		y[n2-j] = y1 + wi*y2;
		sum += wr*y2;
	};
	realft(y, m, 1);
	y[2] = sum;
  	for (j = 4; j <= n; j += 2) {
      		sum += y[j]; y[j] = sum;
	};
	if (isign  ==  -1) {
		even = y[1]; odd = y[2];
      		for (j = 3; j <= n-1; j += 2) {
	  		even += y[j]; odd += y[j+1];
		};
      		enf0 = 2.0*(even - odd);
      		sumo = y[1] - enf0;
      		sume = (2.0*odd/n) - sumo;
      		y[1] = 0.5*enf0;
      		y[2] -= sume;
      		for (j = 3; j <= n-1; j += 2) {
	  		y[j] -= sumo;
	  		y[j+1] -= sume;
		};
    	};
}
//--------------------------------------------------------------------