src/framework/LinearInterpolation.h

0001 /*******************************************************************************
0002  * Copyright (c) The JETSCAPE Collaboration, 2018
0003  *
0004  * Modular, task-based framework for simulating all aspects of heavy-ion collisions
0005  *
0006  * For the list of contributors see AUTHORS.
0007  *
0008  * Report issues at https://github.com/JETSCAPE/JETSCAPE/issues
0009  *
0010  * or via email to bugs.jetscape@gmail.com
0011  *
0012  * Distributed under the GNU General Public License 3.0 (GPLv3 or later).
0013  * See COPYING for details.
0014  ******************************************************************************/
0015
0016 #ifndef LINEARINTERPOLATION_H
0017 #define LINEARINTERPOLATION_H
0018
0019 #include <cmath>
0020
0021 #include "RealType.h"
0022
0023 namespace Jetscape {
0024
0025 /// any type with + and scale * overloaded can use this function
0026 template <class type>
0027 type LinearInt(real x0, real x1, type y0, type y1, real x) {
0028   type temp = ((x - x0) * y1 + (x1 - x) * y0) / (x1 - x0);
0029   return temp;
0030 }
0031
0032 // inspired by numerical recipes
0033 // x0,x1: grid points in x-direction
0034 // y0,y1: grid points in y-direction
0035 // f0-f3: function value starting at x0,y0, continue counterclockwise
0036 // put differently: f0=f(x0,y0)
0037 // f1=f(x1,y0)
0038 // f2=f(x1,y1)
0039 // f3=f(x0,y1)
0040 template <class type>
0041 type BilinearInt(real x0, real x1, real y0, real y1, type f0, type f1, type f2,
0042                  type f3, real x, real y) {
0043   type temp;
0044   real t = (x - x0) / (x1 - x0);
0045   real u = (y - y0) / (y1 - y0);
0046   if ((std::isfinite(u) == 1) && (std::isfinite(t) == 1)) {
0047     temp = (1 - t) * (1 - u) * f0 + t * (1 - u) * f1 + t * u * f2 +
0048            (1 - t) * u * f3;
0049   } else {
0050     if (std::isfinite(u) == 0)
0051       temp = LinearInt(x0, x1, f0, f2, x);
0052     if (std::isfinite(t) == 0)
0053       temp = LinearInt(y0, y1, f0, f2, y);
0054   }
0055   return temp;
0056 }
0057
0058 // 3D linear interpolation
0059 template <class type>
0060 type TrilinearInt(real x0, real x1, real y0, real y1, real z0, real z1,
0061                   type f000, type f001, type f010, type f011, type f100,
0062                   type f101, type f110, type f111, real x, real y, real z) {
0063   type temp;
0064   real t = (x - x0) / (x1 - x0);
0065   real u = (y - y0) / (y1 - y0);
0066   real v = (z - z0) / (z1 - z0);
0067
0068   if ((std::isfinite(u) == 1) && (std::isfinite(t) == 1) &&
0069       (std::isfinite(v) == 1)) {
0070     temp = (1 - t) * (1 - u) * (1 - v) * f000;
0071     temp = temp + (1 - t) * (1 - u) * v * f001;
0072     temp = temp + (1 - t) * u * (1 - v) * f010;
0073     temp = temp + (1 - t) * u * v * f011;
0074     temp = temp + t * (1 - u) * (1 - v) * f100;
0075     temp = temp + t * (1 - u) * v * f101;
0076     temp = temp + t * u * (1 - v) * f110;
0077     temp = temp + t * u * v * f111;
0078   } else {
0079     if (std::isfinite(t) == 0)
0080       temp = BilinearInt(y0, y1, z0, z1, f000, f010, f011, f001, y, z);
0081     if (std::isfinite(v) == 0)
0082       temp = BilinearInt(x0, x1, y0, y1, f000, f100, f110, f010, x, y);
0083     if (std::isfinite(u) == 0)
0084       temp = BilinearInt(x0, x1, z0, z1, f000, f100, f101, f001, x, z);
0085   }
0086   return temp;
0087 }
0088
0089 } //end namespace Jetscape
0090
0091 #endif // LINEARINTERPOLATION_H