src/tools/correl.c

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