]> git.mxchange.org Git - flightgear.git/blob - src/WeatherCM/sphrintp.cpp
The code to find the highest hit below you didn't work quite right when
[flightgear.git] / 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 //---------------------------------------------------------------------------