/** N-dimensional fast Fourier transform using zero-based arrays.
      * Cooley-Tukey FFT algorithm; see Press et al Numerical Recipes 1986.
      * data[0,1,2...] = {real0, imag0, real1, imag1, ....} in and out.
      * ndim is how many dimensions, = 1 for 1-dimensional transform.
      * nn[idim] is how many complex points in each dimension.
      * Each nn must be an exact power of two.
      *
      * @author M.Lampton Java edition (c) 2002 Space Sciences Lab UC Berkeley
      */
    static public void fourn(double data[], int nn[], int ndim, int isign)
    {
        int idim;
        int i1, i2, i3, i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
        int ibit, k1, k2, n, nprev, nrem, ntot;
        double tempi, tempr;
        double theta, wi, wpi, wpr, wr, wtemp;

        for (ntot=1, idim=0; idim<ndim; idim++)
          ntot *= nn[idim]; 
        nprev=1;
        for (idim=ndim-1; idim>=0; idim--)
        {
            n = nn[idim]; 
            nrem = ntot/(n*nprev);
            ip1 = nprev*2;  
            ip2 = ip1*n;  
            ip3 = ip2*nrem;
            i2rev = 0;   
            for (i2=0; 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;  
                            double t = data[i3];    
                            data[i3] = data[i3rev]; 
                            data[i3rev] = t;      
                            t = data[i3+1];   
                            data[i3+1] = data[i3rev+1]; 
                            data[i3rev+1] = t;         
                        }
                    }
                }
                ibit = ip2/2;  
                while (ibit>=ip1 && i2rev>=ibit)    
                {
                    i2rev -= ibit;
                    ibit /= 2; 
                }
                i2rev += ibit;
            }
            ifp1=ip1;
            while (ifp1 < ip2)
            {
                ifp2 = ifp1*2; 
                theta = isign*2*3.14159265358979324/(ifp2/ip1);
                wtemp = Math.sin(0.5*theta);
                wpr = -2.0*wtemp*wtemp;
                wpi = Math.sin(theta);
                wr = 1.0;
                wi=0.0;
                for (i3=0; 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;
        }
        if (isign < 0)
          for (i1=0; i1<2*ntot; i1++)
            data[i1] /= ntot;       
    }