src/gmxlib/fourn.c

   1 /*
   2  * $Id$
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROningen Mixture of Alchemy and Childrens' Stories
  35  */
  36 /* This file is completely threadsafe - keep it that way! */
  37 #include <math.h>
  38 #include "typedefs.h"
  39 #include "nr.h"
  40
  41 #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
  42
  43 void fourn(real data[],int nn[],int ndim,int isign)
  44 {
  45   int i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2;
  46   int ibit,idim,k1,k2,n,nprev,nrem,ntot;
  47   real tempi,tempr;
  48   double theta,wi,wpi,wpr,wr,wtemp;
  49
  50   ntot=1;
  51   for (idim=1;idim<=ndim;idim++)
  52     ntot *= nn[idim];
  53   nprev=1;
  54   for (idim=ndim;idim>=1;idim--) {
  55     n=nn[idim];
  56     nrem=ntot/(n*nprev);
  57     ip1=nprev << 1;
  58     ip2=ip1*n;
  59     ip3=ip2*nrem;
  60     i2rev=1;
  61     for (i2=1;i2<=ip2;i2+=ip1) {
  62       if (i2 < i2rev) {
  63         for (i1=i2;i1<=i2+ip1-2;i1+=2) {
  64           for (i3=i1;i3<=ip3;i3+=ip2) {
  65             i3rev=i2rev+i3-i2;
  66             SWAP(data[i3],data[i3rev]);
  67             SWAP(data[i3+1],data[i3rev+1]);
  68           }
  69         }
  70       }
  71       ibit=ip2 >> 1;
  72       while (ibit >= ip1 && i2rev > ibit) {
  73         i2rev -= ibit;
  74         ibit >>= 1;
  75       }
  76       i2rev += ibit;
  77     }
  78     ifp1=ip1;
  79     while (ifp1 < ip2) {
  80       ifp2=ifp1 << 1;
  81       theta=isign*6.28318530717959/(ifp2/ip1);
  82       wtemp=sin(0.5*theta);
  83       wpr = -2.0*wtemp*wtemp;
  84       wpi=sin(theta);
  85       wr=1.0;
  86       wi=0.0;
  87       for (i3=1;i3<=ifp1;i3+=ip1) {
  88         for (i1=i3;i1<=i3+ip1-2;i1+=2) {
  89           for (i2=i1;i2<=ip3;i2+=ifp2) {
  90             k1=i2;
  91             k2=k1+ifp1;
  92             tempr=wr*data[k2]-wi*data[k2+1];
  93             tempi=wr*data[k2+1]+wi*data[k2];
  94             data[k2]=data[k1]-tempr;
  95             data[k2+1]=data[k1+1]-tempi;
  96             data[k1] += tempr;
  97             data[k1+1] += tempi;
  98           }
  99         }
 100         wr=(wtemp=wr)*wpr-wi*wpi+wr;
 101         wi=wi*wpr+wtemp*wpi+wi;
 102       }
 103       ifp1=ifp2;
 104     }
 105     nprev *= n;
 106   }
 107 }
 108
 109 #undef SWAP