src/WeatherCM/sphrintp.cpp

   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 <simgear/compiler.h>
  33 #include STL_IOSTREAM
  34 #include <simgear/debug/logstream.hxx>
  35
  36 #include <math.h>
  37 #include "sphrintp.h"
  38
  39 SG_USING_NAMESPACE(std);
  40 SG_USING_STD(cout);
  41
  42
  43 static const double PI = 4.0*atan(1.0);
  44 static const double TWOPI = 2.0*PI;
  45
  46 //---------------------------------------------------------------------------
  47 SphereInterpolate::SphereInterpolate (int n, const double* x,
  48                                          const double* y, const double* z,
  49                                          const unsigned int* f)
  50 {
  51     // Assumes (x[i],y[i],z[i]) is unit length for all 0 <= i < n.
  52     // For complete spherical coverage, include the two antipodal points
  53     // (0,0,1,f(0,0,1)) and (0,0,-1,f(0,0,-1)) in the data set.
  54
  55     SG_LOG(SG_MATH, SG_DEBUG, "Initialising spherical interpolator.");
  56     SG_LOG(SG_MATH, SG_DEBUG, "[  0%] Allocating memory");
  57
  58     theta = new double[3*n];
  59     phi = new double[3*n];
  60     func = new unsigned int[3*n];
  61
  62     // convert data to spherical coordinates
  63     int i;
  64
  65     for (i = 0; i < n; i++)
  66     {
  67         GetSphericalCoords(x[i],y[i],z[i],theta[i],phi[i]);
  68         func[i] = f[i];
  69     }
  70
  71     // use periodicity to get wrap-around in the Delaunay triangulation
  72     SG_LOG(SG_MATH, SG_DEBUG, "[ 10%] copying vertices for wrap-aroundr");
  73     int j, k;
  74     for (i = 0, j = n, k = 2*n; i < n; i++, j++, k++)
  75     {
  76         theta[j] = theta[i]+TWOPI;
  77         theta[k] = theta[i]-TWOPI;
  78         phi[j] = phi[i];
  79         phi[k] = phi[i];
  80         func[j] = func[i];
  81         func[k] = func[i];
  82     }
  83
  84     pInterp = new mgcLinInterp2D(3*n,theta,phi,func);
  85
  86     SG_LOG(SG_MATH, SG_DEBUG, "[100%] Finished initialising spherical interpolator.");
  87 }
  88
  89 SphereInterpolate::SphereInterpolate (int n, const sgVec2* p, const unsigned int* f)
  90 {
  91     // Assumes (x[i],y[i],z[i]) is unit length for all 0 <= i < n.
  92     // For complete spherical coverage, include the two antipodal points
  93     // (0,0,1,f(0,0,1)) and (0,0,-1,f(0,0,-1)) in the data set.
  94     SG_LOG(SG_MATH, SG_DEBUG, "Initialising spherical interpolator.");
  95     SG_LOG(SG_MATH, SG_DEBUG, "[  0%] Allocating memory");
  96
  97     theta = new double[3*n];
  98     phi = new double[3*n];
  99     func = new unsigned int[3*n];
 100
 101     // convert data to spherical coordinates
 102     SG_LOG(SG_MATH, SG_DEBUG, "[ 10%] copying vertices for wrap-around");
 103
 104     int i, j, k;
 105     for (i = 0, j = n, k = 2*n; i < n; i++, j++, k++)
 106     {
 107         phi[i] = p[i][0];
 108         theta[i] = p[i][1];
 109         func[i] = f[i];
 110
 111         // use periodicity to get wrap-around in the Delaunay triangulation
 112         phi[j] = phi[i];
 113         phi[k] = phi[i];
 114         theta[j] = theta[i]+TWOPI;
 115         theta[k] = theta[i]-TWOPI;
 116         func[j] = func[i];
 117         func[k] = func[i];
 118     }
 119
 120     pInterp = new mgcLinInterp2D(3*n,theta,phi,func);
 121
 122     SG_LOG(SG_MATH, SG_DEBUG, "[100%] Finished initialising spherical interpolator.");
 123 }
 124 //---------------------------------------------------------------------------
 125 SphereInterpolate::~SphereInterpolate ()
 126 {
 127     delete pInterp;
 128     delete[] theta;
 129     delete[] phi;
 130     delete[] func;
 131 }
 132 //---------------------------------------------------------------------------
 133 void SphereInterpolate::GetSphericalCoords (const double x, const double y, const double z,
 134                                                double& thetaAngle,
 135                                                double& phiAngle) const
 136 {
 137     // Assumes (x,y,z) is unit length.  Returns -PI <= thetaAngle <= PI
 138     // and 0 <= phiAngle <= PI.
 139
 140     if ( z < 1.0f )
 141     {
 142         if ( z > -1.0f )
 143         {
 144             thetaAngle = atan2(y,x);
 145             phiAngle = acos(z);
 146         }
 147         else
 148         {
 149             thetaAngle = -PI;
 150             phiAngle = PI;
 151         }
 152     }
 153     else
 154     {
 155         thetaAngle = -PI;
 156         phiAngle = 0.0f;
 157     }
 158 }
 159 //---------------------------------------------------------------------------
 160 int SphereInterpolate::Evaluate (const double x, const double y, const double z, EvaluateData& f) const
 161 {
 162     // assumes (x,y,z) is unit length
 163
 164     double thetaAngle, phiAngle;
 165     GetSphericalCoords(x,y,z,thetaAngle,phiAngle);
 166     return pInterp->Evaluate(thetaAngle,phiAngle,f);
 167 }
 168 //---------------------------------------------------------------------------
 169 int SphereInterpolate::Evaluate (const double thetaAngle, const double phiAngle, EvaluateData& f) const
 170 {
 171     return pInterp->Evaluate(thetaAngle,phiAngle,f);
 172 }
 173 //---------------------------------------------------------------------------