src/tools/correl.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2008, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  42
  43 #include <stdio.h>
  44 #include <stdlib.h>
  45 #include <math.h>
  46 #include "gmx_fft.h"
  47 #include "smalloc.h"
  48 #include "correl.h"
  49
  50 #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
  51
  52 void four1(real data[],int nn,int isign)
  53 {
  54   int n,mmax,m,j,istep,i;
  55   double wtemp,wr,wpr,wpi,wi,theta;
  56   real tempr,tempi;
  57
  58   n=nn << 1;
  59   j=1;
  60   for (i=1;i<n;i+=2) {
  61     if (j > i) {
  62       SWAP(data[j],data[i]);
  63       SWAP(data[j+1],data[i+1]);
  64     }
  65     m=n >> 1;
  66     while (m >= 2 && j > m) {
  67       j -= m;
  68       m >>= 1;
  69     }
  70     j += m;
  71   }
  72   mmax=2;
  73   while (n > mmax) {
  74     istep=2*mmax;
  75     theta=6.28318530717959/(isign*mmax);
  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 (m=1;m<mmax;m+=2) {
  82       for (i=m;i<=n;i+=istep) {
  83         j=i+mmax;
  84         tempr=wr*data[j]-wi*data[j+1];
  85         tempi=wr*data[j+1]+wi*data[j];
  86         data[j]=data[i]-tempr;
  87         data[j+1]=data[i+1]-tempi;
  88         data[i] += tempr;
  89         data[i+1] += tempi;
  90       }
  91       wr=(wtemp=wr)*wpr-wi*wpi+wr;
  92       wi=wi*wpr+wtemp*wpi+wi;
  93     }
  94     mmax=istep;
  95   }
  96 }
  97
  98 #undef SWAP
  99
 100 static void realft(real data[],int n,int isign)
 101 {
 102   int i,i1,i2,i3,i4,n2p3;
 103   real c1=0.5,c2,h1r,h1i,h2r,h2i;
 104   double wr,wi,wpr,wpi,wtemp,theta;
 105
 106   theta=3.141592653589793/(double) n;
 107   if (isign == 1) {
 108     c2 = -0.5;
 109     four1(data,n,1);
 110   } else {
 111     c2=0.5;
 112     theta = -theta;
 113   }
 114   wtemp=sin(0.5*theta);
 115   wpr = -2.0*wtemp*wtemp;
 116   wpi=sin(theta);
 117   wr=1.0+wpr;
 118   wi=wpi;
 119   n2p3=2*n+3;
 120   for (i=2;i<=n/2;i++) {
 121     i4=1+(i3=n2p3-(i2=1+(i1=i+i-1)));
 122     h1r=c1*(data[i1]+data[i3]);
 123     h1i=c1*(data[i2]-data[i4]);
 124     h2r = -c2*(data[i2]+data[i4]);
 125     h2i=c2*(data[i1]-data[i3]);
 126     data[i1]=h1r+wr*h2r-wi*h2i;
 127     data[i2]=h1i+wr*h2i+wi*h2r;
 128     data[i3]=h1r-wr*h2r+wi*h2i;
 129     data[i4] = -h1i+wr*h2i+wi*h2r;
 130     wr=(wtemp=wr)*wpr-wi*wpi+wr;
 131     wi=wi*wpr+wtemp*wpi+wi;
 132   }
 133   if (isign == 1) {
 134     data[1] = (h1r=data[1])+data[2];
 135     data[2] = h1r-data[2];
 136   } else {
 137     data[1]=c1*((h1r=data[1])+data[2]);
 138     data[2]=c1*(h1r-data[2]);
 139     four1(data,n,-1);
 140   }
 141 }
 142
 143 static void twofft(real data1[],real data2[],real fft1[],real fft2[],int n)
 144 {
 145   int nn3,nn2,jj,j;
 146   real rep,rem,aip,aim;
 147
 148   nn3=1+(nn2=2+n+n);
 149   for (j=1,jj=2;j<=n;j++,jj+=2) {
 150     fft1[jj-1]=data1[j];
 151     fft1[jj]=data2[j];
 152   }
 153   four1(fft1,n,1);
 154   fft2[1]=fft1[2];
 155   fft1[2]=fft2[2]=0.0;
 156   for (j=3;j<=n+1;j+=2) {
 157     rep=0.5*(fft1[j]+fft1[nn2-j]);
 158     rem=0.5*(fft1[j]-fft1[nn2-j]);
 159     aip=0.5*(fft1[j+1]+fft1[nn3-j]);
 160     aim=0.5*(fft1[j+1]-fft1[nn3-j]);
 161     fft1[j]=rep;
 162     fft1[j+1]=aim;
 163     fft1[nn2-j]=rep;
 164     fft1[nn3-j] = -aim;
 165     fft2[j]=aip;
 166     fft2[j+1] = -rem;
 167     fft2[nn2-j]=aip;
 168     fft2[nn3-j]=rem;
 169   }
 170 }
 171
 172 void correl(real data1[],real data2[],int n,real ans[])
 173 {
 174   int     no2,i;
 175   real dum,*fft;
 176
 177   snew(fft,2*n+1);
 178   twofft(data1,data2,fft,ans,n);
 179   no2=n/2;
 180   for (i=2;i<=n+2;i+=2) {
 181     dum      = ans[i-1];
 182     ans[i-1] = (fft[i-1]*dum+fft[i]*ans[i])/no2;
 183     ans[i]   = (fft[i]*dum-fft[i-1]*ans[i])/no2;
 184   }
 185   ans[2]=ans[n+1];
 186   realft(ans,no2,-1);
 187   sfree(fft);
 188 }