| 44 | 
  | 
 * | 
| 45 | 
  | 
 *  Created by Xiuquan Sun on 05/09/06. | 
| 46 | 
  | 
 *  @author  Xiuquan Sun  | 
| 47 | 
< | 
 *  @version $Id: Hxy.cpp,v 1.2 2006-05-16 02:06:37 gezelter Exp $ | 
| 47 | 
> | 
 *  @version $Id: Hxy.cpp,v 1.3 2006-05-16 20:38:23 gezelter Exp $ | 
| 48 | 
  | 
 * | 
| 49 | 
  | 
 */ | 
| 50 | 
  | 
 | 
| 60 | 
  | 
#include<string.h> | 
| 61 | 
  | 
#include<stdlib.h> | 
| 62 | 
  | 
#include<math.h> | 
| 63 | 
– | 
#ifndef WITHOUT_FFTW | 
| 64 | 
– | 
#include<fftw3.h> | 
| 65 | 
– | 
#endif | 
| 63 | 
  | 
 | 
| 64 | 
  | 
namespace oopse { | 
| 65 | 
  | 
   | 
| 87 | 
  | 
  } | 
| 88 | 
  | 
 | 
| 89 | 
  | 
  void Hxy::process() { | 
| 90 | 
< | 
#ifndef WITHOUT_FFTW   | 
| 94 | 
< | 
 | 
| 90 | 
> | 
#if defined(HAVE_FFTW_H) || defined(HAVE_DFFTW_H) || defined(HAVE_FFTW3_H) | 
| 91 | 
  | 
    DumpReader reader(info_, dumpFilename_);     | 
| 92 | 
  | 
    int nFrames = reader.getNFrames(); | 
| 93 | 
  | 
    nProcessed_ = nFrames/step_; | 
| 111 | 
  | 
      reader.readFrame(istep); | 
| 112 | 
  | 
      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot(); | 
| 113 | 
  | 
      nMolecules = info_->getNGlobalMolecules(); | 
| 114 | 
< | 
 | 
| 114 | 
> | 
       | 
| 115 | 
  | 
      Mat3x3d hmat = currentSnapshot_->getHmat(); | 
| 120 | 
– | 
 | 
| 116 | 
  | 
       | 
| 117 | 
< | 
      fftw_complex *in, *out; | 
| 117 | 
> | 
#ifdef HAVE_FFTW3_H | 
| 118 | 
  | 
      fftw_plan p; | 
| 119 | 
+ | 
#else | 
| 120 | 
+ | 
      fftwnd_plan p; | 
| 121 | 
+ | 
#endif | 
| 122 | 
+ | 
      fftw_complex *in, *out; | 
| 123 | 
+ | 
 | 
| 124 | 
  | 
       | 
| 125 | 
  | 
      in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (nBinsX_*nBinsY_)); | 
| 126 | 
  | 
      out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *(nBinsX_*nBinsY_)); | 
| 127 | 
< | 
      p =  fftw_plan_dft_2d(nBinsX_,  | 
| 128 | 
< | 
                            nBinsY_,  | 
| 129 | 
< | 
                            in, out, | 
| 130 | 
< | 
                            FFTW_FORWARD,  | 
| 131 | 
< | 
                            FFTW_ESTIMATE); | 
| 127 | 
> | 
 | 
| 128 | 
> | 
#ifdef HAVE_FFTW3_H | 
| 129 | 
> | 
      p = fftw_plan_dft_2d(nBinsX_, nBinsY_, in, out,  | 
| 130 | 
> | 
                           FFTW_FORWARD, FFTW_ESTIMATE);  | 
| 131 | 
> | 
#else | 
| 132 | 
> | 
      p = fftw2d_create_plan(nBinsX_, nBinsY_, FFTW_FORWARD, FFTW_ESTIMATE); | 
| 133 | 
> | 
#endif | 
| 134 | 
  | 
       | 
| 135 | 
  | 
      int i, j;    | 
| 136 | 
  | 
       | 
| 267 | 
  | 
          c_im(in[newindex]) = 0.0; | 
| 268 | 
  | 
        }  | 
| 269 | 
  | 
      } | 
| 270 | 
< | 
       | 
| 271 | 
< | 
      fftw_execute(p);  | 
| 270 | 
> | 
#ifdef HAVE_FFTW3_H | 
| 271 | 
> | 
      fftw_execute(p); | 
| 272 | 
> | 
#else | 
| 273 | 
> | 
      fftwnd_one(p, in, out); | 
| 274 | 
> | 
#endif | 
| 275 | 
  | 
       | 
| 276 | 
  | 
      for (i=0; i< nBinsX_; i++) { | 
| 277 | 
  | 
        for(j=0; j< nBinsY_; j++) { | 
| 279 | 
  | 
          mag[newindex] = pow(c_re(out[newindex]),2) + pow(c_im(out[newindex]),2); | 
| 280 | 
  | 
        } | 
| 281 | 
  | 
      } | 
| 282 | 
< | 
       | 
| 282 | 
> | 
#ifdef HAVE_FFTW3_H | 
| 283 | 
  | 
      fftw_destroy_plan(p); | 
| 284 | 
+ | 
#else | 
| 285 | 
+ | 
      fftwnd_destroy_plan(p); | 
| 286 | 
+ | 
#endif       | 
| 287 | 
  | 
      fftw_free(out); | 
| 288 | 
  | 
      fftw_free(in); | 
| 289 | 
  | 
 |