fourier.C

#include "random.h"
#include <math.h>

/*
    GenRan: Non-Gaussian Random field generator.
    Copyright (C) 2006 Peter Grassl
    Version 1.0-Devel

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
*/

//Definitions for four1
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

void four1(double data[], int nn, int isign)
  /*Replaces data[1..2*nn] by its discrete Fourier transform if isign is input as 1; or replaces data [1..2*nn] by nn times its inverse discrete Fourier transform, if isign is input as -1. data is a complex array of length nn or, equivalently, a real array of length 2*nn. nn MUST be an integer power of 2 (this is not checked for!). */
{
  unsigned long n, mmax,m,j,istep,i;
  double wtemp,wr,wpr,wpi,wi,theta;
  double tempr, tempi;
  
  n=nn << 1;
  j=1;
  for (i=1;i<n;i+=2) {
    if (j > i){
      SWAP(data[j],data[i]);
      SWAP(data[j+1],data[i+1]);
    }
    m =nn;
    while(m >= 2 && j > m) {
      j -= m;
      m >>= 1;
    }
    j +=m;
  }
  //Here begins the Danielson-Lanczos section of the routine
  mmax=2;
  while(n > mmax){
    istep=mmax << 1;
    theta=isign*(6.28318530717959/mmax);
    wtemp=sin(0.5*theta);
    wpr = -2.0*wtemp*wtemp;
    wpi=sin(theta);
    wr=1.0;
    wi=0.0;
    for(m=1;m<mmax;m+=2){
      for(i=m;i<=n;i+=istep){
        j=i+mmax;
        tempr=wr*data[j]-wi*data[j+1];
        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;
  }
}


// void  fourn(float data[], unsigned long nn[], int ndim, int isign)
// /*Replaces data by itndim-dimensional discrete Fourier transform, if isign is input as 1. nn[1..ndim] is an integer array containing the lengths of each dimension (number of complex values), which MUST all be powers of 2. data is a real array of length twice the product of these lengths, in which the data are stored as in a multidimensional complex array: real and imaginary parts of each element are in consecutive locations, and the rightmost index of the array increases most rapidly as one proceeds along the data. For a two-dimensional, this is equivalent to storing the array by rows. If isign is input as -1, data its inverse transform times the product of the lengths of all dimensions. 
// */
// {
//   int idim;
//   unsigned long i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2;
//   unsigned long ibit,k1,k2,n,nprev,nrem,ntot;
//   float tempi,tempr;
//   double theta,wi,wpi,wpr,wr,wtemp;
//   for (ntot=1,idim=1;idim<=1;idim++)
//     ntot *= nn[idim];
//   nprev=1;
  
//   for (ntot=1,idim=1;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=(float)wr*data[k2]-(float)wi*data[k2+1];
//             tempi=(float)wr*data[k2+1]+(float)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;
//   }
// }