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