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