src/snark/foru_frfilt.cpp

   1 /*
   2  ***********************************************************
   3  $SNARK_Header: S N A R K  1 4 - A PICTURE RECONSTRUCTION PROGRAM $
   4  $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/foru_frfilt.cpp $
   5  $LastChangedRevision: 85 $
   6  $Date: 2014-07-02 16:07:08 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9  */
  10
  11 #include <cstdio>
  12 #include <cmath>
  13
  14 #include "blkdta.h"
  15 #include "geom.h"
  16 #include "consts.h"
  17 #include "sinc.h"
  18
  19 #include "foru.h"
  20
  21 void foru_class::frfilt(REAL* vector, INTEGER size1, INTEGER nfiltr, REAL cutoff, REAL phi)
  22 {
  23         REAL wvqntm;
  24         REAL wvnum;
  25         INTEGER k;
  26         REAL fact;
  27
  28         wvqntm = (REAL) 1.0 / ((REAL) (size1 / 2));
  29         if (GeoPar.vri)
  30                 wvqntm /= MAX0((REAL) fabs(cos(phi)), (REAL) fabs(sin(phi)));
  31         wvnum = 0.;
  32         switch (nfiltr)
  33         {
  34         case 0:
  35                 for (k = 0; k <= size1 - 2; k += 2)
  36                 {  // changed "<" to "<=" to match Fortran code. Lajos, Jan 13, 2005
  37                         if (wvnum > cutoff)
  38                         {
  39                                 vector[k] = 0.;
  40                                 vector[k + 1] = 0.;
  41                         }
  42                         wvnum += wvqntm;
  43                 }
  44                 return;
  45
  46         case 1:
  47                 for (k = 0; k <= size1 - 2; k += 2)
  48                 {  // changed "<" to "<=" to match Fortran code. Lajos, Jan 13, 2005
  49                         if (wvnum > cutoff)
  50                         {
  51                                 vector[k] = 0.;
  52                                 vector[k + 1] = 0.;
  53                         }
  54                         else
  55                         {
  56                                 fact = sinc(wvnum / cutoff);
  57                                 vector[k] *= fact;
  58                                 vector[k + 1] *= fact;
  59                         }
  60                         wvnum += wvqntm;
  61                 }
  62                 return;
  63
  64         case 2:
  65                 for (k = 0; k <= size1 - 2; k += 2)
  66                 {  // changed "<" to "<=" to match Fortran code. Lajos, Jan 13, 2005
  67                         if (wvnum > cutoff)
  68                         {
  69                                 vector[k] = 0.;
  70                                 vector[k + 1] = 0.;
  71                         }
  72                         else
  73                         {
  74                                 fact = (REAL) cos(Consts.pid2 * wvnum / cutoff);
  75                                 vector[k] *= fact;
  76                                 vector[k + 1] *= fact;
  77                         }
  78                         wvnum += wvqntm;
  79                 }
  80                 return;
  81
  82         case 3:
  83                 for (k = 0; k <= size1 - 2; k += 2)
  84                 {  // changed "<" to "<=" to match Fortran code. Lajos, Jan 13, 2005
  85                         if (wvnum > 1.)
  86                         {
  87                                 vector[k] = 0.;
  88                                 vector[k + 1] = 0.;
  89                         }
  90                         else
  91                         {
  92                                 fact = cutoff
  93                                                 + ((REAL) 1.0 - cutoff) * (REAL) cos(Consts.pi * wvnum);
  94                                 vector[k] *= fact;
  95                                 vector[k + 1] *= fact;
  96                         }
  97                         wvnum += wvqntm;
  98                 }
  99                 return;
 100         }
 101 }