Merge pull request #5 from atrettin/avg_ramp

Averaging with linear ramp smoothing

include/SQuIDS/SUNalg.h

@@ -530,14 +530,45 @@ class SU_vector{
   ///               greater than zero allows for a gradual transition to cut-off, rather
   ///               than applying a hard step function.
   void LowPassFilter(double* buffer, double cutoff, double scale) const{
+    if (std::abs(scale) > std::abs(cutoff))
+      throw std::runtime_error("Linear ramp scale cannot be larger than cutoff.");
     auto& suv1=*this;
     size_t offset=GetEvolveBufferSize()/2;
     double* CX=buffer;
     double* SX=buffer+offset;
-    double freq;
+    double term;
     int i;
 #include "SU_inc/LowPassFilterSelect.txt"    
   }
+
+  ///\brief Compute averaging filter for pre-computed sine and cosine evaluations.
+  ///
+  /// This function applies an averaging filter to the pre-evolution buffer computed by
+  /// PrepareEvolve to filter out high frequencies in the same manner as the averaging
+  /// mechanism, meaning that it applies a cut-off to the number of rotations rather
+  /// than frequency. The difference to the already existing mechanism is that this
+  /// filter can be applied as a post-processing step (like the low-pass filter), which
+  /// means that it can also be applied after averaging over a distance. This filter is
+  /// NOT appropriate for the RHS of the state density evolution! Use the low-pass
+  /// filter instead for that purpose.
+  ///\param buffer  The buffer with evaluated sine/cosine values from PrepareEvolve.
+  ///\param t       Evolution time.
+  ///\param cutoff  Cut-off frequency of the filter. Sine and cosine evaluations 
+  ///               with input (frequencies * time) higher than this will be set to zero
+  ///\param scale   Distance in (frequency * time) between cut-off and pass-through. A
+  ///               value of greater than zero allows for a gradual transition to
+  ///               cut-off, rather than applying a hard step function.
+  void AvgRampFilter(double* buffer, double t, double cutoff, double scale) const{
+    if (std::abs(scale) > std::abs(cutoff))
+      throw std::runtime_error("Linear ramp scale cannot be larger than cutoff.");
+    auto& suv1=*this;
+    size_t offset=GetEvolveBufferSize()/2;
+    double* CX=buffer;
+    double* SX=buffer+offset;
+    double term;
+    int i;
+#include "SU_inc/AvgRampFilterSelect.txt"
+  }
 
 
   ///\brief Compute the time evolution of the SU_vector

include/SQuIDS/SU_inc/ApplyLowPassFilter.txt → include/SQuIDS/SU_inc/ApplyLowPassRamp.txt

@@ -23,13 +23,13 @@
  *         atrettin@icecube.wisc.edu                                           *
  ******************************************************************************/
  // greater than cutoff implies setting to zero
- if (fabs(freq) > fabs(cutoff)){
+ if (fabs(term) > fabs(cutoff)){
     SX[i] = 0;
     CX[i] = 0;
  // If we are below the cutoff, but above the threshold where we begin to apply the 
  // filter gradually, we do that.
- }else if(fabs(freq) > fabs(cutoff) - fabs(scale)){
-    SX[i] *= (fabs(cutoff) - fabs(freq))/fabs(scale);
-    CX[i] *= (fabs(cutoff) - fabs(freq))/fabs(scale);
+ }else if(fabs(term) > fabs(cutoff) - fabs(scale)){
+    SX[i] *= (fabs(cutoff) - fabs(term))/fabs(scale);
+    CX[i] *= (fabs(cutoff) - fabs(term))/fabs(scale);
  }
  // if neither of the cases above applied, we do nothing

include/SQuIDS/SU_inc/AvgRampFilterSelect.txt

@@ -0,0 +1,43 @@
+ /******************************************************************************
+ *    This program is free software: you can redistribute it and/or modify     *
+ *   it under the terms of the GNU General Public License as published by      *
+ *   the Free Software Foundation, either version 3 of the License, or         *
+ *   (at your option) any later version.                                       *
+ *                                                                             *
+ *   This program is distributed in the hope that it will be useful,           *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of            *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             *
+ *   GNU General Public License for more details.                              *
+ *                                                                             *
+ *   You should have received a copy of the GNU General Public License         *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.     *
+ *                                                                             *   
+ *   Authors:                                                                  *
+ *      Carlos Arguelles (University of Wisconsin Madison)                     * 
+ *         carguelles@icecube.wisc.edu                                         *
+ *      Jordi Salvado (University of Wisconsin Madison)                        *
+ *         jsalvado@icecube.wisc.edu                                           *
+ *      Christopher Weaver (University of Wisconsin Madison)                   * 
+ *         cweaver@icecube.wisc.edu                                            *
+ *      Alexander Trettin (DESY)                                               *
+ *         atrettin@icecube.wisc.edu                                           *
+ ******************************************************************************/
+switch (dim){
+  case 2: 
+    #include "AvgWithRampSU2.txt"
+    break;
+  case 3: 
+    #include "AvgWithRampSU3.txt"
+    break;
+  case 4: 
+    #include "AvgWithRampSU4.txt"
+    break;
+  case 5: 
+    #include "AvgWithRampSU5.txt"
+    break;
+  case 6: 
+    #include "AvgWithRampSU6.txt"
+    break;
+  default:
+    throw std::runtime_error("SU_vector Evolution Error: Only dimensions up to " SQUIDS_MAX_HILBERT_DIM_STR " are supported");
+}

include/SQuIDS/SU_inc/AvgWithRampSU2.txt

@@ -0,0 +1,27 @@
+ /******************************************************************************
+ *    This program is free software: you can redistribute it and/or modify     *
+ *   it under the terms of the GNU General Public License as published by      *
+ *   the Free Software Foundation, either version 3 of the License, or         *
+ *   (at your option) any later version.                                       *
+ *                                                                             *
+ *   This program is distributed in the hope that it will be useful,           *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of            *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             *
+ *   GNU General Public License for more details.                              *
+ *                                                                             *
+ *   You should have received a copy of the GNU General Public License         *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.     *
+ *                                                                             *   
+ *   Authors:                                                                  *
+ *      Carlos Arguelles (University of Wisconsin Madison)                     * 
+ *         carguelles@icecube.wisc.edu                                         *
+ *      Jordi Salvado (University of Wisconsin Madison)                        *
+ *         jsalvado@icecube.wisc.edu                                           *
+ *      Christopher Weaver (University of Wisconsin Madison)                   * 
+ *         cweaver@icecube.wisc.edu                                            *
+ *      Alexander Trettin (DESY)                                               *
+ *         atrettin@icecube.wisc.edu                                           *
+ ******************************************************************************/
+term=2*t*suv1.components[3];
+i = 0;
+#include "ApplyLowPassRamp.txt"

include/SQuIDS/SU_inc/AvgWithRampSU3.txt

@@ -0,0 +1,38 @@
+ /******************************************************************************
+ *    This program is free software: you can redistribute it and/or modify     *
+ *   it under the terms of the GNU General Public License as published by      *
+ *   the Free Software Foundation, either version 3 of the License, or         *
+ *   (at your option) any later version.                                       *
+ *                                                                             *
+ *   This program is distributed in the hope that it will be useful,           *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of            *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             *
+ *   GNU General Public License for more details.                              *
+ *                                                                             *
+ *   You should have received a copy of the GNU General Public License         *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.     *
+ *                                                                             *   
+ *   Authors:                                                                  *
+ *      Carlos Arguelles (University of Wisconsin Madison)                     * 
+ *         carguelles@icecube.wisc.edu                                         *
+ *      Jordi Salvado (University of Wisconsin Madison)                        *
+ *         jsalvado@icecube.wisc.edu                                           *
+ *      Christopher Weaver (University of Wisconsin Madison)                   * 
+ *         cweaver@icecube.wisc.edu                                            *
+ *      Alexander Trettin (DESY)                                               *
+ *         atrettin@icecube.wisc.edu                                           *
+ ******************************************************************************/
+term=2*t*suv1.components[4];
+i = 0;
+#include "ApplyLowPassRamp.txt"
+//printf("0: term: %4.2E   scale:  %4.2E  time: %4.2E \n",term,scale,t);
+
+term=t*(suv1.components[4] + sqrt(3)*suv1.components[8]);
+i = 1;
+#include "ApplyLowPassRamp.txt"
+//printf("1: term: %4.2E   scale:  %4.2E  time: %4.2E \n",term,scale,t);
+
+term=t*(suv1.components[4] - sqrt(3)*suv1.components[8]);
+i = 2;
+#include "ApplyLowPassRamp.txt"
+//printf("2: term: %4.2E   scale:  %4.2E  time: %4.2E \n",term,scale,t);

include/SQuIDS/SU_inc/AvgWithRampSU4.txt

@@ -0,0 +1,47 @@
+ /******************************************************************************
+ *    This program is free software: you can redistribute it and/or modify     *
+ *   it under the terms of the GNU General Public License as published by      *
+ *   the Free Software Foundation, either version 3 of the License, or         *
+ *   (at your option) any later version.                                       *
+ *                                                                             *
+ *   This program is distributed in the hope that it will be useful,           *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of            *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             *
+ *   GNU General Public License for more details.                              *
+ *                                                                             *
+ *   You should have received a copy of the GNU General Public License         *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.     *
+ *                                                                             *   
+ *   Authors:                                                                  *
+ *      Carlos Arguelles (University of Wisconsin Madison)                     * 
+ *         carguelles@icecube.wisc.edu                                         *
+ *      Jordi Salvado (University of Wisconsin Madison)                        *
+ *         jsalvado@icecube.wisc.edu                                           *
+ *      Christopher Weaver (University of Wisconsin Madison)                   * 
+ *         cweaver@icecube.wisc.edu                                            *
+ *      Alexander Trettin (DESY)                                               *
+ *         atrettin@icecube.wisc.edu                                           *
+ ******************************************************************************/
+term=2*t*suv1.components[5];
+i = 0;
+#include "ApplyLowPassRamp.txt"
+
+term=t*(suv1.components[5] + sqrt(3)*suv1.components[10]);
+i = 1;
+#include "ApplyLowPassRamp.txt"
+
+term=t*suv1.components[5] + (t*(suv1.components[10] + 2*sqrt(2)*suv1.components[15]))/sqrt(3);
+i = 2;
+#include "ApplyLowPassRamp.txt"
+
+term=t*(suv1.components[5] - sqrt(3)*suv1.components[10]);
+i = 3;
+#include "ApplyLowPassRamp.txt"
+
+term=t*suv1.components[5] - (t*(suv1.components[10] + 2*sqrt(2)*suv1.components[15]))/sqrt(3);
+i = 4;
+#include "ApplyLowPassRamp.txt"
+
+term=(2*t*(suv1.components[10] - sqrt(2)*suv1.components[15]))/sqrt(3);
+i = 5;
+#include "ApplyLowPassRamp.txt"

include/SQuIDS/SU_inc/AvgWithRampSU5.txt

@@ -0,0 +1,63 @@
+ /******************************************************************************
+ *    This program is free software: you can redistribute it and/or modify     *
+ *   it under the terms of the GNU General Public License as published by      *
+ *   the Free Software Foundation, either version 3 of the License, or         *
+ *   (at your option) any later version.                                       *
+ *                                                                             *
+ *   This program is distributed in the hope that it will be useful,           *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of            *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             *
+ *   GNU General Public License for more details.                              *
+ *                                                                             *
+ *   You should have received a copy of the GNU General Public License         *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.     *
+ *                                                                             *   
+ *   Authors:                                                                  *
+ *      Carlos Arguelles (University of Wisconsin Madison)                     * 
+ *         carguelles@icecube.wisc.edu                                         *
+ *      Jordi Salvado (University of Wisconsin Madison)                        *
+ *         jsalvado@icecube.wisc.edu                                           *
+ *      Christopher Weaver (University of Wisconsin Madison)                   * 
+ *         cweaver@icecube.wisc.edu                                            *
+ *      Alexander Trettin (DESY)                                               *
+ *         atrettin@icecube.wisc.edu                                           *
+ ******************************************************************************/
+term=2*t*suv1.components[6];
+i = 0;
+#include "ApplyLowPassRamp.txt"
+
+term=t*(suv1.components[6] + sqrt(3)*suv1.components[12]);
+i = 1;
+#include "ApplyLowPassRamp.txt"
+
+term=t*suv1.components[6] + (t*(suv1.components[12] + 2*sqrt(2)*suv1.components[18]))/sqrt(3);
+i = 2;
+#include "ApplyLowPassRamp.txt"
+
+term=t*suv1.components[6] + (t*suv1.components[12])/sqrt(3) + (t*suv1.components[18])/sqrt(6) + sqrt(2.5)*t*suv1.components[24];
+i = 3;
+#include "ApplyLowPassRamp.txt"
+
+term=t*(suv1.components[6] - sqrt(3)*suv1.components[12]);
+i = 4;
+#include "ApplyLowPassRamp.txt"
+
+term=t*suv1.components[6] - (t*(suv1.components[12] + 2*sqrt(2)*suv1.components[18]))/sqrt(3);
+i = 5;
+#include "ApplyLowPassRamp.txt"
+
+term=t*suv1.components[6] - (t*(2*sqrt(3)*suv1.components[12] + sqrt(6)*suv1.components[18] + 3*sqrt(10)*suv1.components[24]))/6.;
+i = 6;
+#include "ApplyLowPassRamp.txt"
+
+term=(2*t*(suv1.components[12] - sqrt(2)*suv1.components[18]))/sqrt(3);
+i = 7;
+#include "ApplyLowPassRamp.txt"
+
+term=(t*(4*sqrt(3)*suv1.components[12] - sqrt(6)*suv1.components[18] - 3*sqrt(10)*suv1.components[24]))/6.;
+i = 8;
+#include "ApplyLowPassRamp.txt"
+
+term=sqrt(1.5)*t*suv1.components[18] - sqrt(2.5)*t*suv1.components[24];
+i = 9;
+#include "ApplyLowPassRamp.txt"

include/SQuIDS/SU_inc/AvgWithRampSU6.txt

@@ -0,0 +1,84 @@
+ /******************************************************************************
+ *    This program is free software: you can redistribute it and/or modify     *
+ *   it under the terms of the GNU General Public License as published by      *
+ *   the Free Software Foundation, either version 3 of the License, or         *
+ *   (at your option) any later version.                                       *
+ *                                                                             *
+ *   This program is distributed in the hope that it will be useful,           *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of            *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             *
+ *   GNU General Public License for more details.                              *
+ *                                                                             *
+ *   You should have received a copy of the GNU General Public License         *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.     *
+ *                                                                             *   
+ *   Authors:                                                                  *
+ *      Carlos Arguelles (University of Wisconsin Madison)                     * 
+ *         carguelles@icecube.wisc.edu                                         *
+ *      Jordi Salvado (University of Wisconsin Madison)                        *
+ *         jsalvado@icecube.wisc.edu                                           *
+ *      Christopher Weaver (University of Wisconsin Madison)                   * 
+ *         cweaver@icecube.wisc.edu                                            *
+ *      Alexander Trettin (DESY)                                               *
+ *         atrettin@icecube.wisc.edu                                           *
+ ******************************************************************************/
+
+term=2*t*suv1.components[7];
+i = 0;
+#include "ApplyLowPassRamp.txt"
+
+term=t*(suv1.components[7] + sqrt(3)*suv1.components[14]);
+i = 1;
+#include "ApplyLowPassRamp.txt"
+
+term=t*suv1.components[7] + (t*(suv1.components[14] + 2*sqrt(2)*suv1.components[21]))/sqrt(3);
+i = 2;
+#include "ApplyLowPassRamp.txt"
+
+term=t*suv1.components[7] + (t*suv1.components[14])/sqrt(3) + (t*suv1.components[21])/sqrt(6) + sqrt(2.5)*t*suv1.components[28];
+i = 3;
+#include "ApplyLowPassRamp.txt"
+
+term=t*suv1.components[7] + (t*suv1.components[14])/sqrt(3) + (t*suv1.components[21])/sqrt(6) + (t*suv1.components[28])/sqrt(10) + 2*sqrt(0.6)*t*suv1.components[35];
+i = 4;
+#include "ApplyLowPassRamp.txt"
+
+term=t*(suv1.components[7] - sqrt(3)*suv1.components[14]);
+i = 5;
+#include "ApplyLowPassRamp.txt"
+
+term=t*suv1.components[7] - (t*(suv1.components[14] + 2*sqrt(2)*suv1.components[21]))/sqrt(3);
+i = 6;
+#include "ApplyLowPassRamp.txt"
+
+term=t*suv1.components[7] - (t*(2*sqrt(3)*suv1.components[14] + sqrt(6)*suv1.components[21] + 3*sqrt(10)*suv1.components[28]))/6.;
+i = 7;
+#include "ApplyLowPassRamp.txt"
+
+term=t*suv1.components[7] - (t*(10*sqrt(3)*suv1.components[14] + 5*sqrt(6)*suv1.components[21] + 3*sqrt(10)*suv1.components[28] + 12*sqrt(15)*suv1.components[35]))/30.;
+i = 8;
+#include "ApplyLowPassRamp.txt"
+
+term=(2*t*(suv1.components[14] - sqrt(2)*suv1.components[21]))/sqrt(3);
+i = 9;
+#include "ApplyLowPassRamp.txt"
+
+term=(t*(4*sqrt(3)*suv1.components[14] - sqrt(6)*suv1.components[21] - 3*sqrt(10)*suv1.components[28]))/6.;
+i = 10;
+#include "ApplyLowPassRamp.txt"
+
+term=(2*t*suv1.components[14])/sqrt(3) - (t*suv1.components[21])/sqrt(6) - (t*suv1.components[28])/sqrt(10) - 2*sqrt(0.6)*t*suv1.components[35];
+i = 11;
+#include "ApplyLowPassRamp.txt"
+
+term=sqrt(1.5)*t*suv1.components[21] - sqrt(2.5)*t*suv1.components[28];
+i = 12;
+#include "ApplyLowPassRamp.txt"
+
+term=sqrt(1.5)*t*suv1.components[21] - (t*suv1.components[28])/sqrt(10) - 2*sqrt(0.6)*t*suv1.components[35];
+i = 13;
+#include "ApplyLowPassRamp.txt"
+
+term=2*sqrt(0.4)*t*suv1.components[28] - 2*sqrt(0.6)*t*suv1.components[35];
+i = 14;
+#include "ApplyLowPassRamp.txt"