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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Green Red Orange Magenta Azure Cyan Skyblue
  28  */
  29 static char *SRCID_fourn_c = "$Id$";
  30
  31 #include <math.h>
  32 #include "typedefs.h"
  33 #include "nr.h"
  34
  35 #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
  36
  37 void fourn(real data[],int nn[],int ndim,int isign)
  38 {
  39   int i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2;
  40   int ibit,idim,k1,k2,n,nprev,nrem,ntot;
  41   real tempi,tempr;
  42   double theta,wi,wpi,wpr,wr,wtemp;
  43
  44   ntot=1;
  45   for (idim=1;idim<=ndim;idim++)
  46     ntot *= nn[idim];
  47   nprev=1;
  48   for (idim=ndim;idim>=1;idim--) {
  49     n=nn[idim];
  50     nrem=ntot/(n*nprev);
  51     ip1=nprev << 1;
  52     ip2=ip1*n;
  53     ip3=ip2*nrem;
  54     i2rev=1;
  55     for (i2=1;i2<=ip2;i2+=ip1) {
  56       if (i2 < i2rev) {
  57         for (i1=i2;i1<=i2+ip1-2;i1+=2) {
  58           for (i3=i1;i3<=ip3;i3+=ip2) {
  59             i3rev=i2rev+i3-i2;
  60             SWAP(data[i3],data[i3rev]);
  61             SWAP(data[i3+1],data[i3rev+1]);
  62           }
  63         }
  64       }
  65       ibit=ip2 >> 1;
  66       while (ibit >= ip1 && i2rev > ibit) {
  67         i2rev -= ibit;
  68         ibit >>= 1;
  69       }
  70       i2rev += ibit;
  71     }
  72     ifp1=ip1;
  73     while (ifp1 < ip2) {
  74       ifp2=ifp1 << 1;
  75       theta=isign*6.28318530717959/(ifp2/ip1);
  76       wtemp=sin(0.5*theta);
  77       wpr = -2.0*wtemp*wtemp;
  78       wpi=sin(theta);
  79       wr=1.0;
  80       wi=0.0;
  81       for (i3=1;i3<=ifp1;i3+=ip1) {
  82         for (i1=i3;i1<=i3+ip1-2;i1+=2) {
  83           for (i2=i1;i2<=ip3;i2+=ifp2) {
  84             k1=i2;
  85             k2=k1+ifp1;
  86             tempr=wr*data[k2]-wi*data[k2+1];
  87             tempi=wr*data[k2+1]+wi*data[k2];
  88             data[k2]=data[k1]-tempr;
  89             data[k2+1]=data[k1+1]-tempi;
  90             data[k1] += tempr;
  91             data[k1+1] += tempi;
  92           }
  93         }
  94         wr=(wtemp=wr)*wpr-wi*wpi+wr;
  95         wi=wi*wpr+wtemp*wpi+wi;
  96       }
  97       ifp1=ifp2;
  98     }
  99     nprev *= n;
 100   }
 101 }
 102
 103 #undef SWAP