2 WARNING - Do not remove this header.
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.
9 Dave Eberly says that the conditions for use are:
11 * You may distribute the original source code to others at no charge.
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.
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
23 * The original code comes with absolutely no warranty and no guarantee is
24 made that the code is bug-free.
26 This does not seem incompatible with GPL - so this modified version
27 is hereby placed under GPL along with the rest of FlightGear.
32 #include <simgear/compiler.h>
38 SG_USING_NAMESPACE(std);
42 static const double PI = 4.0*atan(1.0);
43 static const double TWOPI = 2.0*PI;
45 //---------------------------------------------------------------------------
46 SphereInterpolate::SphereInterpolate (int n, const double* x,
47 const double* y, const double* z,
48 const unsigned int* f)
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.
54 cout << "Initialising spherical interpolator.\n";
55 cout << "[ 0%] Allocating memory \r";
57 theta = new double[3*n];
58 phi = new double[3*n];
59 func = new unsigned int[3*n];
61 // convert data to spherical coordinates
64 for (i = 0; i < n; i++)
66 GetSphericalCoords(x[i],y[i],z[i],theta[i],phi[i]);
70 // use periodicity to get wrap-around in the Delaunay triangulation
71 cout << "[ 10%] copying vertices for wrap-around\r";
73 for (i = 0, j = n, k = 2*n; i < n; i++, j++, k++)
75 theta[j] = theta[i]+TWOPI;
76 theta[k] = theta[i]-TWOPI;
83 pInterp = new mgcLinInterp2D(3*n,theta,phi,func);
85 cout << "[100%] Finished initialising spherical interpolator. \n";
88 SphereInterpolate::SphereInterpolate (int n, const sgVec2* p, const unsigned int* f)
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";
96 theta = new double[3*n];
97 phi = new double[3*n];
98 func = new unsigned int[3*n];
100 // convert data to spherical coordinates
101 cout << "[ 10%] copying vertices for wrap-around \r";
104 for (i = 0, j = n, k = 2*n; i < n; i++, j++, k++)
110 // use periodicity to get wrap-around in the Delaunay triangulation
113 theta[j] = theta[i]+TWOPI;
114 theta[k] = theta[i]-TWOPI;
119 pInterp = new mgcLinInterp2D(3*n,theta,phi,func);
121 cout << "[100%] Finished initialising spherical interpolator. \n";
123 //---------------------------------------------------------------------------
124 SphereInterpolate::~SphereInterpolate ()
131 //---------------------------------------------------------------------------
132 void SphereInterpolate::GetSphericalCoords (const double x, const double y, const double z,
134 double& phiAngle) const
136 // Assumes (x,y,z) is unit length. Returns -PI <= thetaAngle <= PI
137 // and 0 <= phiAngle <= PI.
143 thetaAngle = atan2(y,x);
158 //---------------------------------------------------------------------------
159 int SphereInterpolate::Evaluate (const double x, const double y, const double z, EvaluateData& f) const
161 // assumes (x,y,z) is unit length
163 double thetaAngle, phiAngle;
164 GetSphericalCoords(x,y,z,thetaAngle,phiAngle);
165 return pInterp->Evaluate(thetaAngle,phiAngle,f);
167 //---------------------------------------------------------------------------
168 int SphereInterpolate::Evaluate (const double thetaAngle, const double phiAngle, EvaluateData& f) const
170 return pInterp->Evaluate(thetaAngle,phiAngle,f);
172 //---------------------------------------------------------------------------