1 files changed, 82 insertions, 0 deletions
diff --git a/externals/grill/vasp/source/rdx2fft.cpp b/externals/grill/vasp/source/rdx2fft.cpp
new file mode 100644
index 00000000..b4ecf6c9
--- /dev/null
+++ b/externals/grill/vasp/source/rdx2fft.cpp
@@ -0,0 +1,82 @@
+#include <math.h>
+
+#define PI 3.1415926535897932384f
+
+//////////////////////////////////////////////////////////////////////////
+
+/* calculate bidirectional fourier transform of complex data radix 2 */
+/* adapted from subroutine FOUREA listed in                          */
+/* Programs for Digital Signal Processing                            */
+/* edited by Digital Signal Processing Committee                     */
+/* IEEE Acoustics Speech and Signal Processing Committee             */
+/* Chapter 1 Section 1.1 Page 1.1-4,5                                */
+/* direct -1 forward +1 reverse                                      */
+
+bool fft_bidir_complex_radix2(int size,float *real,float *imag,int direct)
+{
+  int i,j,m,mmax,istep;
+  float c,s,treal,timag,theta; 
+
+  /* compute transform */
+
+  j=1;
+  for(i=1;i<=size;i++)
+  {
+    if(i<j)
+    {
+      treal=real[j-1];
+      timag=imag[j-1];
+      real[j-1]=real[i-1];
+      imag[j-1]=imag[i-1];
+      real[i-1]=treal;
+      imag[i-1]=timag;
+    }
+    m=size/2;
+    while(j>m)
+    {
+      j-=m;
+      m=(m+1)/2;
+    }
+    j+=m;
+  }
+  mmax=1;
+  while(size>mmax)
+  {
+    istep=2*mmax;
+    for(m=1;m<=mmax;m++)
+    {
+      theta=PI*(float)direct*(float)(m-1)/(float)mmax;
+      c=(float)cos(theta);
+      s=(float)sin(theta);
+      for(i=m;i<=size;i+=istep)
+      {
+	j=i+mmax;
+	treal=real[j-1]*c-imag[j-1]*s;
+	timag=imag[j-1]*c+real[j-1]*s;
+	real[j-1]=real[i-1]-treal;
+	imag[j-1]=imag[i-1]-timag;
+	real[i-1]+=treal;
+	imag[i-1]+=timag;
+      }
+    }
+    mmax=istep;
+  }
+
+  return true;
+}
+
+#if 0
+/* calculate forward fourier transform of complex data radix 2 */
+
+bool fft_fwd_complex_radix2(int size,float *real,float *imag)
+{
+  return fft_bidir_complex_radix2(size,real,imag,-1);
+}
+
+/* calculate inverse fourier transform of complex data radix 2 */
+
+bool fft_inv_complex_radix2(int size,float *real,float *imag)
+{
+  return fft_bidir_complex_radix2(size,real,imag,1);
+}
+#endif