ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/staticProps/Hxy.cpp
(Generate patch)

Comparing trunk/src/applications/staticProps/Hxy.cpp (file contents):
Revision 956 by gezelter, Tue May 16 02:06:37 2006 UTC vs.
Revision 963 by tim, Wed May 17 21:51:42 2006 UTC

# Line 44 | Line 44
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.4 2006-05-17 21:51:42 tim Exp $
48   *
49   */
50  
# Line 60 | Line 60
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    
# Line 90 | Line 87 | namespace oopse {
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_;
94      
95 <    std::vector<double> mag, newmag;
96 <    double lenX_, lenY_;
97 <    double gridX_, gridY_;
98 <    double halfBoxX_, halfBoxY_;
95 >    std::vector<RealType> mag, newmag;
96 >    RealType lenX_, lenY_;
97 >    RealType gridX_, gridY_;
98 >    RealType halfBoxX_, halfBoxY_;
99      int binNoX, binNoY;
100 <    double interpsum, value;
100 >    RealType interpsum, value;
101      int ninterp, px, py, newp;
102      int newx, newy, newindex, index;
103      int new_i, new_j, new_index;
104 <    double freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
105 <    double maxfreqx, maxfreqy, maxfreq, dfreq;
104 >    RealType freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
105 >    RealType maxfreqx, maxfreqy, maxfreq, dfreq;
106      int whichbin;
107      int nMolecules;
108      
# Line 115 | Line 111 | namespace oopse {
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        
# Line 156 | Line 158 | namespace oopse {
158        gridX_ = lenX_ /(nBinsX_);
159        gridY_ = lenY_ /(nBinsY_);
160        
161 <      double halfBoxX_ = lenX_ / 2.0;      
162 <      double halfBoxY_ = lenY_ / 2.0;      
161 >      RealType halfBoxX_ = lenX_ / 2.0;      
162 >      RealType halfBoxY_ = lenY_ / 2.0;      
163        
164        if (evaluator_.isDynamic()) {
165          seleMan_.setSelectionSet(evaluator_.evaluate());
# Line 196 | Line 198 | namespace oopse {
198          for(j = 0; j < nBinsY_; j++){
199            newindex = i * nBinsY_ + j;
200            if(gridsample_[newindex] > 0){
201 <            gridZ_[newindex] = gridZ_[newindex] / (double)gridsample_[newindex];
201 >            gridZ_[newindex] = gridZ_[newindex] / (RealType)gridsample_[newindex];
202            }
203          }
204        }
# Line 250 | Line 252 | namespace oopse {
252                ninterp++;
253              }
254          
255 <            value = interpsum / (double)ninterp;
255 >            value = interpsum / (RealType)ninterp;
256              
257              gridZ_[newindex] = value;
258            }
# Line 265 | Line 267 | namespace oopse {
267            c_im(in[newindex]) = 0.0;
268          }
269        }
270 + #ifdef HAVE_FFTW3_H
271 +      fftw_execute(p);
272 + #else
273 +      fftwnd_one(p, in, out);
274 + #endif
275        
269 –      fftw_execute(p);
270 –      
276        for (i=0; i< nBinsX_; i++) {
277          for(j=0; j< nBinsY_; j++) {
278            newindex = i*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  
# Line 325 | Line 333 | namespace oopse {
333        //  printf("%lf\t%lf\t%lf\t%lf\n", dx, dy, maxfreqx, maxfreqy);
334        
335        maxfreq = sqrt(maxfreqx*maxfreqx + maxfreqy*maxfreqy);
336 <      dfreq = maxfreq/(double)(nbins_-1);
336 >      dfreq = maxfreq/(RealType)(nbins_-1);
337      
338        //printf("%lf\n", dfreq);
339        
# Line 335 | Line 343 | namespace oopse {
343        for (i=0; i< nBinsX_; i++) {
344          for(j=0; j< nBinsY_; j++) {
345            
346 <          freq_x = (double)(i - zero_freq_x)*maxfreqx*2 / nBinsX_;
347 <          freq_y = (double)(j - zero_freq_y)*maxfreqy*2 / nBinsY_;
346 >          freq_x = (RealType)(i - zero_freq_x)*maxfreqx*2 / nBinsX_;
347 >          freq_y = (RealType)(j - zero_freq_y)*maxfreqy*2 / nBinsY_;
348            
349            freq = sqrt(freq_x*freq_x + freq_y*freq_y);
350            
# Line 350 | Line 358 | namespace oopse {
358        
359        for ( i = 0; i < nbins_; i++) {
360          if ( samples[i][istep] > 0) {
361 <          bin[i][istep] = 4.0 * sqrt(bin[i][istep] / (double)samples[i][istep]) / (double)nMolecules;
361 >          bin[i][istep] = 4.0 * sqrt(bin[i][istep] / (RealType)samples[i][istep]) / (RealType)nMolecules;
362          }
363        }    
364        
# Line 361 | Line 369 | namespace oopse {
369          sum_bin[i] += bin[i][j];
370          sum_bin_sq[i] += bin[i][j] * bin[i][j];
371        }
372 <      avg_bin[i] = sum_bin[i] / (double)nFrames;
373 <      avg_bin_sq[i] = sum_bin_sq[i] / (double)nFrames;
372 >      avg_bin[i] = sum_bin[i] / (RealType)nFrames;
373 >      avg_bin_sq[i] = sum_bin_sq[i] / (RealType)nFrames;
374        for (int j = 0; j < nFrames; j++) {
375          errbin_sum[i] += pow((bin[i][j] - avg_bin[i]), 2);
376          errbin_sum_sq[i] += pow((bin[i][j] * bin[i][j] - avg_bin_sq[i]), 2);
377        }
378 <      errbin[i] = sqrt( errbin_sum[i] / (double)nFrames );
379 <      errbin_sq[i] = sqrt( errbin_sum_sq[i] / (double)nFrames );
378 >      errbin[i] = sqrt( errbin_sum[i] / (RealType)nFrames );
379 >      errbin_sq[i] = sqrt( errbin_sum_sq[i] / (RealType)nFrames );
380      }
381      
382      printSpectrum();

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines