src/WeatherCM/sphrintp.inl

   1 /*
   2   WARNING - Do not remove this header.
   3
   4   This code is a templated version of the 'magic-software' spherical
   5   interpolation code by Dave Eberly. The original (un-hacked) code can be
   6   obtained from here: http://www.magic-software.com/gr_appr.htm
   7   This code is derived from linintp2.h/cpp and sphrintp.h/cpp.
   8
   9   Dave Eberly says that the conditions for use are:
  10
  11   * You may distribute the original source code to others at no charge.
  12
  13   * You may modify the original source code and distribute it to others at
  14     no charge. The modified code must be documented to indicate that it is
  15     not part of the original package.
  16
  17   * You may use this code for non-commercial purposes. You may also
  18     incorporate this code into commercial packages. However, you may not
  19     sell any of your source code which contains my original and/or modified
  20     source code. In such a case, you need to factor out my code and freely
  21     distribute it.
  22
  23   * The original code comes with absolutely no warranty and no guarantee is
  24     made that the code is bug-free.
  25
  26   This does not seem incompatible with GPL - so this modified version
  27   is hereby placed under GPL along with the rest of FlightGear.
  28
  29                               Christian Mayer
  30 */
  31
  32 #include <math.h>
  33 #include "sphrintp.h"
  34
  35 static const double PI = 4.0*atan(1.0);
  36 static const double TWOPI = 2.0*PI;
  37
  38 //---------------------------------------------------------------------------
  39 template<class T>
  40 SphereInterpolate<T>::SphereInterpolate (int n, const double* x,
  41                                          const double* y, const double* z,
  42                                          const T* f)
  43 {
  44     // Assumes (x[i],y[i],z[i]) is unit length for all 0 <= i < n.
  45     // For complete spherical coverage, include the two antipodal points
  46     // (0,0,1,f(0,0,1)) and (0,0,-1,f(0,0,-1)) in the data set.
  47
  48     cout << "Initialising spherical interpolator.\n";
  49     cout << "[  0%] Allocating memory                                           \r";
  50
  51     theta = new double[3*n];
  52     phi = new double[3*n];
  53     func = new T[3*n];
  54
  55     // convert data to spherical coordinates
  56     int i;
  57     T empty;
  58
  59     for (i = 0; i < n; i++)
  60     {
  61         GetSphericalCoords(x[i],y[i],z[i],theta[i],phi[i]);
  62         func[i] = f[i];
  63     }
  64
  65     // use periodicity to get wrap-around in the Delaunay triangulation
  66     cout << "[ 10%] copying vertices for wrap-around\r";
  67     int j, k;
  68     for (i = 0, j = n, k = 2*n; i < n; i++, j++, k++)
  69     {
  70         theta[j] = theta[i]+TWOPI;
  71         theta[k] = theta[i]-TWOPI;
  72         phi[j] = phi[i];
  73         phi[k] = phi[i];
  74         func[j] = func[i];
  75         func[k] = func[i];
  76     }
  77
  78     pInterp = new mgcLinInterp2D<T>(3*n,theta,phi,func);
  79
  80     cout << "[100%] Finished initialising spherical interpolator.                   \n";
  81 }
  82
  83 template<class T>
  84 SphereInterpolate<T>::SphereInterpolate (int n, const sgVec2* p, const T* f)
  85 {
  86     // Assumes (x[i],y[i],z[i]) is unit length for all 0 <= i < n.
  87     // For complete spherical coverage, include the two antipodal points
  88     // (0,0,1,f(0,0,1)) and (0,0,-1,f(0,0,-1)) in the data set.
  89     cout << "Initialising spherical interpolator.\n";
  90     cout << "[  0%] Allocating memory                                           \r";
  91
  92     theta = new double[3*n];
  93     phi = new double[3*n];
  94     func = new T[3*n];
  95
  96     // convert data to spherical coordinates
  97     cout << "[ 10%] copying vertices for wrap-around                              \r";
  98
  99     int i, j, k;
 100     for (i = 0, j = n, k = 2*n; i < n; i++, j++, k++)
 101     {
 102         phi[i] = p[i][0];
 103         theta[i] = p[i][1];
 104         func[i] = f[i];
 105
 106         // use periodicity to get wrap-around in the Delaunay triangulation
 107         phi[j] = phi[i];
 108         phi[k] = phi[i];
 109         theta[j] = theta[i]+TWOPI;
 110         theta[k] = theta[i]-TWOPI;
 111         func[j] = func[i];
 112         func[k] = func[i];
 113     }
 114
 115     pInterp = new mgcLinInterp2D<T>(3*n,theta,phi,func);
 116
 117     cout << "[100%] Finished initialising spherical interpolator.                   \n";
 118 }
 119 //---------------------------------------------------------------------------
 120 template<class T>
 121 SphereInterpolate<T>::~SphereInterpolate ()
 122 {
 123     delete pInterp;
 124     delete[] theta;
 125     delete[] phi;
 126     delete[] func;
 127 }
 128 //---------------------------------------------------------------------------
 129 template<class T>
 130 void SphereInterpolate<T>::GetSphericalCoords (const double x, const double y, const double z,
 131                                                double& thetaAngle,
 132                                                double& phiAngle) const
 133 {
 134     // Assumes (x,y,z) is unit length.  Returns -PI <= thetaAngle <= PI
 135     // and 0 <= phiAngle <= PI.
 136
 137     if ( z < 1.0f )
 138     {
 139         if ( z > -1.0f )
 140         {
 141             thetaAngle = atan2(y,x);
 142             phiAngle = acos(z);
 143         }
 144         else
 145         {
 146             thetaAngle = -PI;
 147             phiAngle = PI;
 148         }
 149     }
 150     else
 151     {
 152         thetaAngle = -PI;
 153         phiAngle = 0.0f;
 154     }
 155 }
 156 //---------------------------------------------------------------------------
 157 template<class T>
 158 int SphereInterpolate<T>::Evaluate (const double x, const double y, const double z, T& f) const
 159 {
 160     // assumes (x,y,z) is unit length
 161
 162     double thetaAngle, phiAngle;
 163     GetSphericalCoords(x,y,z,thetaAngle,phiAngle);
 164     return pInterp->Evaluate(thetaAngle,phiAngle,f);
 165 }
 166 //---------------------------------------------------------------------------
 167 template<class T>
 168 int SphereInterpolate<T>::Evaluate (const double thetaAngle, const double phiAngle, T& f) const
 169 {
 170     return pInterp->Evaluate(thetaAngle,phiAngle,f);
 171 }
 172 //---------------------------------------------------------------------------