--- /dev/null
+/*
+ WARNING - Do not remove this header.
+
+ This code is a templated version of the 'magic-software' spherical
+ interpolation code by Dave Eberly. The original (un-hacked) code can be
+ obtained from here: http://www.magic-software.com/gr_appr.htm
+ This code is derived from linintp2.h/cpp and sphrintp.h/cpp.
+
+ Dave Eberly says that the conditions for use are:
+
+ * You may distribute the original source code to others at no charge.
+
+ * You may modify the original source code and distribute it to others at
+ no charge. The modified code must be documented to indicate that it is
+ not part of the original package.
+
+ * You may use this code for non-commercial purposes. You may also
+ incorporate this code into commercial packages. However, you may not
+ sell any of your source code which contains my original and/or modified
+ source code. In such a case, you need to factor out my code and freely
+ distribute it.
+
+ * The original code comes with absolutely no warranty and no guarantee is
+ made that the code is bug-free.
+
+ This does not seem incompatible with GPL - so this modified version
+ is hereby placed under GPL along with the rest of FlightGear.
+
+ Christian Mayer
+*/
+
+#ifndef LININTP2_H
+#define LININTP2_H
+
+template<class T>
+class mgcLinInterp2D
+{
+public:
+ mgcLinInterp2D (int _numPoints, double* x, double* y, T* _f);
+
+ ~mgcLinInterp2D ();
+
+ double XMin () { return xmin; }
+ double XMax () { return xmax; }
+ double XRange () { return xmax-xmin; }
+ double YMin () { return ymin; }
+ double YMax () { return ymax; }
+ double YRange () { return ymax-ymin; }
+
+ int PointCount () { return numPoints; }
+ void GetPoint (int i, double& x, double& y);
+
+ int EdgeCount () { return numEdges; }
+ void GetEdge (int i, double& x0, double& y0, double& x1, double& y1);
+
+ int TriangleCount () { return numTriangles; }
+ void GetTriangle (int i, double& x0, double& y0, double& x1, double& y1,
+ double& x2, double& y2);
+
+ int Evaluate (double x, double y, T& F);
+
+private:
+ typedef struct
+ {
+ double x, y;
+ }
+ Vertex;
+
+ typedef struct
+ {
+ int vertex[3]; // listed in counterclockwise order
+
+ int adj[3];
+ // adj[0] points to triangle sharing edge (vertex[0],vertex[1])
+ // adj[1] points to triangle sharing edge (vertex[1],vertex[2])
+ // adj[2] points to triangle sharing edge (vertex[2],vertex[0])
+ }
+ Triangle;
+
+ typedef struct
+ {
+ int vertex[2];
+ int triangle[2];
+ int index[2];
+ }
+ Edge;
+
+ int numPoints;
+ double** point;
+ double** tmppoint;
+ T* f;
+
+ double xmin, xmax, ymin, ymax;
+
+
+ int numEdges;
+ Edge* edge;
+
+ int numTriangles;
+ Triangle* triangle;
+
+ int Delaunay2D ();
+ void ComputeBarycenter (Vertex& v0, Vertex& v1, Vertex& v2, Vertex& ver,
+ double c[3]);
+ int InTriangle (Vertex& v0, Vertex& v1, Vertex& v2, Vertex& test);
+};
+
+#include "linintp2.inl"
+
+#endif
--- /dev/null
+/*
+ WARNING - Do not remove this header.
+
+ This code is a templated version of the 'magic-software' spherical
+ interpolation code by Dave Eberly. The original (un-hacked) code can be
+ obtained from here: http://www.magic-software.com/gr_appr.htm
+ This code is derived from linintp2.h/cpp and sphrintp.h/cpp.
+
+ Dave Eberly says that the conditions for use are:
+
+ * You may distribute the original source code to others at no charge.
+
+ * You may modify the original source code and distribute it to others at
+ no charge. The modified code must be documented to indicate that it is
+ not part of the original package.
+
+ * You may use this code for non-commercial purposes. You may also
+ incorporate this code into commercial packages. However, you may not
+ sell any of your source code which contains my original and/or modified
+ source code. In such a case, you need to factor out my code and freely
+ distribute it.
+
+ * The original code comes with absolutely no warranty and no guarantee is
+ made that the code is bug-free.
+
+ This does not seem incompatible with GPL - so this modified version
+ is hereby placed under GPL along with the rest of FlightGear.
+
+ Christian Mayer
+*/
+
+#include <float.h>
+#include <math.h>
+#include <stdlib.h>
+#include "linintp2.h"
+
+//---------------------------------------------------------------------------
+template<class T>
+mgcLinInterp2D<T>::mgcLinInterp2D (int _numPoints, double* x, double* y,
+ T* _f)
+{
+ if ( (numPoints = _numPoints) < 3 )
+ {
+ point = 0;
+ edge = 0;
+ triangle = 0;
+ numTriangles = 0;
+ return;
+ }
+
+ cout << "[ 20%] allocating memory \r";
+
+ point = new double*[numPoints];
+ tmppoint = new double*[numPoints+3];
+ f = new T[numPoints];
+ int i;
+ for (i = 0; i < numPoints; i++)
+ point[i] = new double[2];
+ for (i = 0; i < numPoints+3; i++)
+ tmppoint[i] = new double[2];
+ for (i = 0; i < numPoints; i++)
+ {
+ point[i][0] = tmppoint[i][0] = x[i];
+ point[i][1] = tmppoint[i][1] = y[i];
+
+ f[i] = _f[i];
+ }
+
+ cout << "[ 30%] creating delaunay diagram \r";
+
+ Delaunay2D();
+}
+//---------------------------------------------------------------------------
+template<class T>
+mgcLinInterp2D<T>::~mgcLinInterp2D ()
+{
+ if ( numPoints < 3 )
+ return;
+
+ int i;
+
+ if ( point )
+ {
+ for (i = 0; i < numPoints; i++)
+ delete[] point[i];
+ delete[] point;
+ }
+ if ( tmppoint )
+ {
+ for (i = 0; i < numPoints+3; i++)
+ delete[] tmppoint[i];
+ delete[] tmppoint;
+ }
+
+ delete[] f;
+ delete[] edge;
+ delete[] triangle;
+}
+//---------------------------------------------------------------------------
+template<class T>
+void mgcLinInterp2D<T>::ComputeBarycenter (Vertex& v0, Vertex& v1, Vertex& v2,
+ Vertex& ver, double c[3])
+{
+ double A0 = v0.x-v2.x, B0 = v0.y-v2.y;
+ double A1 = v1.x-v2.x, B1 = v1.y-v2.y;
+ double A2 = ver.x-v2.x, B2 = ver.y-v2.y;
+
+ double m00 = A0*A0+B0*B0, m01 = A0*A1+B0*B1, m11 = A1*A1+B1*B1;
+ double r0 = A2*A0+B2*B0, r1 = A2*A1+B2*B1;
+ double det = m00*m11-m01*m01;
+
+ c[0] = (m11*r0-m01*r1)/det;
+ c[1] = (m00*r1-m01*r0)/det;
+ c[2] = 1-c[0]-c[1];
+}
+//---------------------------------------------------------------------------
+template<class T>
+int mgcLinInterp2D<T>::InTriangle (Vertex& v0, Vertex& v1, Vertex& v2,
+ Vertex& test)
+{
+ const double eps = 1e-08;
+ double tx, ty, nx, ny;
+
+ // test against normal to first edge
+ tx = test.x - v0.x;
+ ty = test.y - v0.y;
+ nx = v0.y - v1.y;
+ ny = v1.x - v0.x;
+ if ( tx*nx + ty*ny < -eps )
+ return 0;
+
+ // test against normal to second edge
+ tx = test.x - v1.x;
+ ty = test.y - v1.y;
+ nx = v1.y - v2.y;
+ ny = v2.x - v1.x;
+ if ( tx*nx + ty*ny < -eps )
+ return 0;
+
+ // test against normal to third edge
+ tx = test.x - v2.x;
+ ty = test.y - v2.y;
+ nx = v2.y - v0.y;
+ ny = v0.x - v2.x;
+ if ( tx*nx + ty*ny < -eps )
+ return 0;
+
+ return 1;
+}
+//---------------------------------------------------------------------------
+template<class T>
+int mgcLinInterp2D<T>::Evaluate (double x, double y, T& F)
+{
+ Vertex ver = { x, y };
+ // determine which triangle contains the target point
+
+ int i;
+ Vertex v0, v1, v2;
+ for (i = 0; i < numTriangles; i++)
+ {
+ Triangle& t = triangle[i];
+ v0.x = point[t.vertex[0]][0];
+ v0.y = point[t.vertex[0]][1];
+ v1.x = point[t.vertex[1]][0];
+ v1.y = point[t.vertex[1]][1];
+ v2.x = point[t.vertex[2]][0];
+ v2.y = point[t.vertex[2]][1];
+
+ if ( InTriangle(v0,v1,v2,ver) )
+ break;
+ }
+
+ if ( i == numTriangles ) // point is outside interpolation region
+ {
+ return 0;
+ }
+
+ Triangle& t = triangle[i]; // (x,y) is in this triangle
+
+ // compute barycentric coordinates with respect to subtriangle
+ double bary[3];
+ ComputeBarycenter(v0,v1,v2,ver,bary);
+
+ // compute barycentric combination of function values at vertices
+ F = bary[0]*f[t.vertex[0]]+bary[1]*f[t.vertex[1]]+bary[2]*f[t.vertex[2]];
+
+ return 1;
+}
+//---------------------------------------------------------------------------
+template<class T>
+int mgcLinInterp2D<T>::Delaunay2D ()
+{
+ int result;
+
+ const double EPSILON = 1e-12;
+ const int TSIZE = 75;
+ const double RANGE = 10.0;
+
+ xmin = tmppoint[0][0];
+ xmax = xmin;
+ ymin = tmppoint[0][1];
+ ymax = ymin;
+
+ int i;
+ for (i = 0; i < numPoints; i++)
+ {
+ double value = tmppoint[i][0];
+ if ( xmax < value )
+ xmax = value;
+ if ( xmin > value )
+ xmin = value;
+
+ value = tmppoint[i][1];
+ if ( ymax < value )
+ ymax = value;
+ if ( ymin > value )
+ ymin = value;
+ }
+
+ double xrange = xmax-xmin, yrange = ymax-ymin;
+ double maxrange = xrange;
+ if ( maxrange < yrange )
+ maxrange = yrange;
+
+ // need to scale the data later to do a correct triangle count
+ double maxrange2 = maxrange*maxrange;
+
+ // tweak the points by very small random numbers
+ double bgs = EPSILON*maxrange;
+ srand(367);
+ for (i = 0; i < numPoints; i++)
+ {
+ tmppoint[i][0] += bgs*(0.5 - rand()/double(RAND_MAX));
+ tmppoint[i][1] += bgs*(0.5 - rand()/double(RAND_MAX));
+ }
+
+ double wrk[2][3] =
+ {
+ { 5*RANGE, -RANGE, -RANGE },
+ { -RANGE, 5*RANGE, -RANGE }
+ };
+ for (i = 0; i < 3; i++)
+ {
+ tmppoint[numPoints+i][0] = xmin+xrange*wrk[0][i];
+ tmppoint[numPoints+i][1] = ymin+yrange*wrk[1][i];
+ }
+
+ int i0, i1, i2, i3, i4, i5, i6, i7, i8, i9, i11;
+ int nts, ii[3];
+ double xx;
+
+ int tsz = 2*TSIZE;
+ int** tmp = new int*[tsz+1];
+ tmp[0] = new int[2*(tsz+1)];
+ for (i0 = 1; i0 < tsz+1; i0++)
+ tmp[i0] = tmp[0] + 2*i0;
+ i1 = 2*(numPoints + 2);
+
+ int* id = new int[i1];
+ for (i0 = 0; i0 < i1; i0++)
+ id[i0] = i0;
+
+ int** a3s = new int*[i1];
+ a3s[0] = new int[3*i1];
+ for (i0 = 1; i0 < i1; i0++)
+ a3s[i0] = a3s[0] + 3*i0;
+ a3s[0][0] = numPoints;
+ a3s[0][1] = numPoints+1;
+ a3s[0][2] = numPoints+2;
+
+ double** ccr = new double*[i1]; // circumscribed centers and radii
+ ccr[0] = new double[3*i1];
+ for (i0 = 1; i0 < i1; i0++)
+ ccr[i0] = ccr[0] + 3*i0;
+ ccr[0][0] = 0.0;
+ ccr[0][1] = 0.0;
+ ccr[0][2] = FLT_MAX;
+
+ nts = 1; // number of triangles
+ i4 = 1;
+
+ cout << "[ 40%] create triangulation \r";
+
+ // compute triangulation
+ for (i0 = 0; i0 < numPoints; i0++)
+ {
+ i1 = i7 = -1;
+ i9 = 0;
+ for (i11 = 0; i11 < nts; i11++)
+ {
+ i1++;
+ while ( a3s[i1][0] < 0 )
+ i1++;
+ xx = ccr[i1][2];
+ for (i2 = 0; i2 < 2; i2++)
+ {
+ double z = tmppoint[i0][i2]-ccr[i1][i2];
+ xx -= z*z;
+ if ( xx < 0 )
+ goto Corner3;
+ }
+ i9--;
+ i4--;
+ id[i4] = i1;
+ for (i2 = 0; i2 < 3; i2++)
+ {
+ ii[0] = 0;
+ if (ii[0] == i2)
+ ii[0]++;
+ for (i3 = 1; i3 < 2; i3++)
+ {
+ ii[i3] = ii[i3-1] + 1;
+ if (ii[i3] == i2)
+ ii[i3]++;
+ }
+ if ( i7 > 1 )
+ {
+ i8 = i7;
+ for (i3 = 0; i3 <= i8; i3++)
+ {
+ for (i5 = 0; i5 < 2; i5++)
+ if ( a3s[i1][ii[i5]] != tmp[i3][i5] )
+ goto Corner1;
+ for (i6 = 0; i6 < 2; i6++)
+ tmp[i3][i6] = tmp[i8][i6];
+ i7--;
+ goto Corner2;
+Corner1:;
+ }
+ }
+ if ( ++i7 > tsz )
+ {
+ // temporary storage exceeded, increase TSIZE
+ result = 0;
+ goto ExitDelaunay;
+ }
+ for (i3 = 0; i3 < 2; i3++)
+ tmp[i7][i3] = a3s[i1][ii[i3]];
+Corner2:;
+ }
+ a3s[i1][0] = -1;
+Corner3:;
+ }
+
+ for (i1 = 0; i1 <= i7; i1++)
+ {
+ for (i2 = 0; i2 < 2; i2++)
+ for (wrk[i2][2] = 0, i3 = 0; i3 < 2; i3++)
+ {
+ wrk[i2][i3] = tmppoint[tmp[i1][i2]][i3]-tmppoint[i0][i3];
+ wrk[i2][2] +=
+ 0.5*wrk[i2][i3]*(tmppoint[tmp[i1][i2]][i3]+
+ tmppoint[i0][i3]);
+ }
+
+ xx = wrk[0][0]*wrk[1][1]-wrk[1][0]*wrk[0][1];
+ ccr[id[i4]][0] = (wrk[0][2]*wrk[1][1]-wrk[1][2]*wrk[0][1])/xx;
+ ccr[id[i4]][1] = (wrk[0][0]*wrk[1][2]-wrk[1][0]*wrk[0][2])/xx;
+
+ for (ccr[id[i4]][2] = 0, i2 = 0; i2 < 2; i2++)
+ {
+ double z = tmppoint[i0][i2]-ccr[id[i4]][i2];
+ ccr[id[i4]][2] += z*z;
+ a3s[id[i4]][i2] = tmp[i1][i2];
+ }
+
+ a3s[id[i4]][2] = i0;
+ i4++;
+ i9++;
+ }
+ nts += i9;
+ }
+
+ // count the number of triangles
+ cout << "[ 50%] count the number of triangles \r";
+
+ numTriangles = 0;
+ i0 = -1;
+ for (i11 = 0; i11 < nts; i11++)
+ {
+ i0++;
+ while ( a3s[i0][0] < 0 )
+ i0++;
+ if ( a3s[i0][0] < numPoints )
+ {
+ for (i1 = 0; i1 < 2; i1++)
+ for (i2 = 0; i2 < 2; i2++)
+ wrk[i1][i2] =
+ tmppoint[a3s[i0][i1]][i2]-tmppoint[a3s[i0][2]][i2];
+
+ if ( fabs(wrk[0][0]*wrk[1][1]-wrk[0][1]*wrk[1][0]) > EPSILON*maxrange2 )
+ numTriangles++;
+ }
+ }
+
+ // create the triangles
+ cout << "[ 60%] create the triangles \r";
+
+ triangle = new Triangle[numTriangles];
+
+ numTriangles = 0;
+ i0 = -1;
+ for (i11 = 0; i11 < nts; i11++)
+ {
+ i0++;
+ while ( a3s[i0][0] < 0 )
+ i0++;
+ if ( a3s[i0][0] < numPoints )
+ {
+ for (i1 = 0; i1 < 2; i1++)
+ for (i2 = 0; i2 < 2; i2++)
+ wrk[i1][i2] =
+ tmppoint[a3s[i0][i1]][i2]-tmppoint[a3s[i0][2]][i2];
+ xx = wrk[0][0]*wrk[1][1]-wrk[0][1]*wrk[1][0];
+ if ( fabs(xx) > EPSILON*maxrange2 )
+ {
+ int delta = xx < 0 ? 1 : 0;
+ Triangle& tri = triangle[numTriangles];
+ tri.vertex[0] = a3s[i0][0];
+ tri.vertex[1] = a3s[i0][1+delta];
+ tri.vertex[2] = a3s[i0][2-delta];
+ tri.adj[0] = -1;
+ tri.adj[1] = -1;
+ tri.adj[2] = -1;
+ numTriangles++;
+ }
+ }
+ }
+
+ // build edge table
+ cout << "[ 70%] build the edge table \r";
+
+ numEdges = 0;
+ edge = new Edge[3*numTriangles];
+
+ int j, j0, j1;
+ for (i = 0; i < numTriangles; i++)
+ {
+ if ( (i%500) == 0)
+ cout << "[ 7" << 10*i/numTriangles << "%] build the edge table \r";
+
+ Triangle& t = triangle[i];
+
+ for (j0 = 0, j1 = 1; j0 < 3; j0++, j1 = (j1+1)%3)
+ {
+ for (j = 0; j < numEdges; j++)
+ {
+ Edge& e = edge[j];
+ if ( (t.vertex[j0] == e.vertex[0]
+ && t.vertex[j1] == e.vertex[1])
+ || (t.vertex[j0] == e.vertex[1]
+ && t.vertex[j1] == e.vertex[0]) )
+ break;
+ }
+ if ( j == numEdges ) // add edge to table
+ {
+ edge[j].vertex[0] = t.vertex[j0];
+ edge[j].vertex[1] = t.vertex[j1];
+ edge[j].triangle[0] = i;
+ edge[j].index[0] = j0;
+ edge[j].triangle[1] = -1;
+ numEdges++;
+ }
+ else // edge already exists, add triangle to table
+ {
+ edge[j].triangle[1] = i;
+ edge[j].index[1] = j0;
+ }
+ }
+ }
+
+ // establish links between adjacent triangles
+ cout << "[ 80%] establishing links between adjacent triangles \r";
+
+ for (i = 0; i < numEdges; i++)
+ {
+ if ( edge[i].triangle[1] != -1 )
+ {
+ j0 = edge[i].triangle[0];
+ j1 = edge[i].triangle[1];
+ triangle[j0].adj[edge[i].index[0]] = j1;
+ triangle[j1].adj[edge[i].index[1]] = j0;
+ }
+ }
+
+ result = 1;
+
+ExitDelaunay:;
+ delete[] tmp[0];
+ delete[] tmp;
+ delete[] id;
+ delete[] a3s[0];
+ delete[] a3s;
+ delete[] ccr[0];
+ delete[] ccr;
+
+ cout << "[ 90%] finsishes delauney triangulation \r";
+
+ return result;
+}
+//---------------------------------------------------------------------------
+template<class T>
+void mgcLinInterp2D<T>::GetPoint (int i, double& x, double& y)
+{
+ // assumes i is valid [can use PointCount() before passing i]
+ x = point[i][0];
+ y = point[i][1];
+}
+//---------------------------------------------------------------------------
+template<class T>
+void mgcLinInterp2D<T>::GetEdge (int i, double& x0, double& y0, double& x1,
+ double& y1)
+{
+ // assumes i is valid [can use EdgeCount() before passing i]
+ int v0 = edge[i].vertex[0], v1 = edge[i].vertex[1];
+
+ x0 = point[v0][0];
+ y0 = point[v0][1];
+ x1 = point[v1][0];
+ y1 = point[v1][1];
+}
+//---------------------------------------------------------------------------
+template<class T>
+void mgcLinInterp2D<T>::GetTriangle (int i, double& x0, double& y0, double& x1,
+ double& y1, double& x2, double& y2)
+{
+ // assumes i is valid [can use TriangleCount() before passing i]
+ int v0 = triangle[i].vertex[0];
+ int v1 = triangle[i].vertex[1];
+ int v2 = triangle[i].vertex[2];
+
+ x0 = point[v0][0];
+ y0 = point[v0][1];
+ x1 = point[v1][0];
+ y1 = point[v1][1];
+ x2 = point[v2][0];
+ y2 = point[v2][1];
+}
+//---------------------------------------------------------------------------
+
--- /dev/null
+/*
+ WARNING - Do not remove this header.
+
+ This code is a templated version of the 'magic-software' spherical
+ interpolation code by Dave Eberly. The original (un-hacked) code can be
+ obtained from here: http://www.magic-software.com/gr_appr.htm
+ This code is derived from linintp2.h/cpp and sphrintp.h/cpp.
+
+ Dave Eberly says that the conditions for use are:
+
+ * You may distribute the original source code to others at no charge.
+
+ * You may modify the original source code and distribute it to others at
+ no charge. The modified code must be documented to indicate that it is
+ not part of the original package.
+
+ * You may use this code for non-commercial purposes. You may also
+ incorporate this code into commercial packages. However, you may not
+ sell any of your source code which contains my original and/or modified
+ source code. In such a case, you need to factor out my code and freely
+ distribute it.
+
+ * The original code comes with absolutely no warranty and no guarantee is
+ made that the code is bug-free.
+
+ This does not seem incompatible with GPL - so this modified version
+ is hereby placed under GPL along with the rest of FlightGear.
+
+ Christian Mayer
+*/
+
+#include <math.h>
+#include "sphrintp.h"
+
+static const double PI = 4.0*atan(1.0);
+static const double TWOPI = 2.0*PI;
+
+//---------------------------------------------------------------------------
+template<class T>
+SphereInterpolate<T>::SphereInterpolate (int n, const double* x,
+ const double* y, const double* z,
+ const T* f)
+{
+ // Assumes (x[i],y[i],z[i]) is unit length for all 0 <= i < n.
+ // For complete spherical coverage, include the two antipodal points
+ // (0,0,1,f(0,0,1)) and (0,0,-1,f(0,0,-1)) in the data set.
+
+ cout << "Initialising spherical interpolator.\n";
+ cout << "[ 0%] Allocating memory \r";
+
+ theta = new double[3*n];
+ phi = new double[3*n];
+ func = new T[3*n];
+
+ // convert data to spherical coordinates
+ int i;
+ T empty;
+
+ for (i = 0; i < n; i++)
+ {
+ GetSphericalCoords(x[i],y[i],z[i],theta[i],phi[i]);
+ func[i] = f[i];
+ }
+
+ // use periodicity to get wrap-around in the Delaunay triangulation
+ cout << "[ 10%] copying vertices for wrap-around\r";
+ int j, k;
+ for (i = 0, j = n, k = 2*n; i < n; i++, j++, k++)
+ {
+ theta[j] = theta[i]+TWOPI;
+ theta[k] = theta[i]-TWOPI;
+ phi[j] = phi[i];
+ phi[k] = phi[i];
+ func[j] = func[i];
+ func[k] = func[i];
+ }
+
+ pInterp = new mgcLinInterp2D<T>(3*n,theta,phi,func);
+
+ cout << "[100%] Finished initialising spherical interpolator. \n";
+}
+
+template<class T>
+SphereInterpolate<T>::SphereInterpolate (int n, const sgVec2* p, const T* f)
+{
+ // Assumes (x[i],y[i],z[i]) is unit length for all 0 <= i < n.
+ // For complete spherical coverage, include the two antipodal points
+ // (0,0,1,f(0,0,1)) and (0,0,-1,f(0,0,-1)) in the data set.
+ cout << "Initialising spherical interpolator.\n";
+ cout << "[ 0%] Allocating memory \r";
+
+ theta = new double[3*n];
+ phi = new double[3*n];
+ func = new T[3*n];
+
+ // convert data to spherical coordinates
+ cout << "[ 10%] copying vertices for wrap-around \r";
+
+ int i, j, k;
+ for (i = 0, j = n, k = 2*n; i < n; i++, j++, k++)
+ {
+ phi[i] = p[i][0];
+ theta[i] = p[i][1];
+ func[i] = f[i];
+
+ // use periodicity to get wrap-around in the Delaunay triangulation
+ phi[j] = phi[i];
+ phi[k] = phi[i];
+ theta[j] = theta[i]+TWOPI;
+ theta[k] = theta[i]-TWOPI;
+ func[j] = func[i];
+ func[k] = func[i];
+ }
+
+ pInterp = new mgcLinInterp2D<T>(3*n,theta,phi,func);
+
+ cout << "[100%] Finished initialising spherical interpolator. \n";
+}
+//---------------------------------------------------------------------------
+template<class T>
+SphereInterpolate<T>::~SphereInterpolate ()
+{
+ delete pInterp;
+ delete[] theta;
+ delete[] phi;
+ delete[] func;
+}
+//---------------------------------------------------------------------------
+template<class T>
+void SphereInterpolate<T>::GetSphericalCoords (const double x, const double y, const double z,
+ double& thetaAngle,
+ double& phiAngle) const
+{
+ // Assumes (x,y,z) is unit length. Returns -PI <= thetaAngle <= PI
+ // and 0 <= phiAngle <= PI.
+
+ if ( z < 1.0f )
+ {
+ if ( z > -1.0f )
+ {
+ thetaAngle = atan2(y,x);
+ phiAngle = acos(z);
+ }
+ else
+ {
+ thetaAngle = -PI;
+ phiAngle = PI;
+ }
+ }
+ else
+ {
+ thetaAngle = -PI;
+ phiAngle = 0.0f;
+ }
+}
+//---------------------------------------------------------------------------
+template<class T>
+int SphereInterpolate<T>::Evaluate (const double x, const double y, const double z, T& f) const
+{
+ // assumes (x,y,z) is unit length
+
+ double thetaAngle, phiAngle;
+ GetSphericalCoords(x,y,z,thetaAngle,phiAngle);
+ return pInterp->Evaluate(thetaAngle,phiAngle,f);
+}
+//---------------------------------------------------------------------------
+template<class T>
+int SphereInterpolate<T>::Evaluate (const double thetaAngle, const double phiAngle, T& f) const
+{
+ return pInterp->Evaluate(thetaAngle,phiAngle,f);
+}
+//---------------------------------------------------------------------------
projects / flightgear.git / commitdiff