CppEphem
test_CEPlanet.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * test_CEPlanet.cpp: CppEphem *
3  * ----------------------------------------------------------------------- *
4  * Copyright © 2019 JCardenzana *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 
22 #include "test_CEPlanet.h"
23 #include "CENamespace.h"
24 
25 
26 /**********************************************************************/
30  CETestSuite()
31 {
32  // Interpolate corrections
34 
35  // Create the date object
37  // Interpret the date object as a julian date
39 
40  // Setup the observer object for observed coordinates
46 }
47 
48 
49 /**********************************************************************/
53 {}
54 
55 
56 /**********************************************************************/
62 {
63  std::cout << "\nTesting CEPlanet:\n";
64 
65  // Run each of the tests
67  test_Planets();
68 
69  return pass();
70 }
71 
72 
73 /**********************************************************************/
79 {
80  CESkyCoord test_coord;
81 
82  // Default constructor
83  CEPlanet test1;
84  test_string(test1.Name(), "undefined", __func__, __LINE__);
85  test(test1.GetCoordinates() == test_coord, __func__, __LINE__);
86 
87  // Default constructors
88  CEPlanet test2("DefaultPlanet", CEAngle::Deg(123.45), CEAngle::Deg(67.89),
90  test_string(test2.Name(), "DefaultPlanet", __func__, __LINE__);
91  test_double(test2.XCoord().Deg(), 123.45, __func__, __LINE__);
92  test_double(test2.YCoord().Deg(), 67.89, __func__, __LINE__);
93  test_int(int(test2.GetCoordSystem()), int(CESkyCoordType::ICRS), __func__, __LINE__);
94 
95  // Copy constructor
96  CEPlanet test3(test2);
97  test_string(test3.Name(), test2.Name(), __func__, __LINE__);
98  test_double(test3.XCoord().Deg(), test2.XCoord().Deg(), __func__, __LINE__);
99  test_double(test3.YCoord().Deg(), test2.YCoord().Deg(), __func__, __LINE__);
100  test_int(int(test3.GetCoordSystem()), int(test2.GetCoordSystem()), __func__, __LINE__);
101 
102  // Copy assignment operator
103  CEPlanet test4 = test2;
104  test_string(test4.Name(), test2.Name(), __func__, __LINE__);
105  test_double(test4.XCoordinate_Deg(), test2.XCoordinate_Deg(), __func__, __LINE__);
106  test_double(test4.YCoordinate_Deg(), test2.YCoordinate_Deg(), __func__, __LINE__);
107  test_int(int(test4.GetCoordSystem()), int(test2.GetCoordSystem()), __func__, __LINE__);
108 
109  return pass();
110 }
111 
112 
113 /**********************************************************************/
118 {
119  // Create the planet object
120  CEPlanet mercury = CEPlanet::Mercury();
121 
122  // Update the coordinates
123  mercury.UpdateCoordinates(base_date_.JD());
124 
125  // Run the actual test with values derived from AstroPy
126  test_planet(mercury,
127  CESkyCoord(CEAngle::Deg(251.1840266663606371),
128  CEAngle::Deg(-25.3023875289787838),
130  CESkyCoord(CEAngle::Deg(197.8036693771968260),
131  CEAngle::Deg(25.7350871794139664),
133  {-0.13721236, -0.4032437, -0.20141521},
134  {0.02137206, -0.00493223, -0.00485005});
135 
136  return pass();
137 }
138 
139 
140 /**********************************************************************/
145 {
146  // Create the planet object
147  CEPlanet venus = CEPlanet::Venus();
148 
149  // Update the coordinates
151 
152  // Run the actual test with values derived from AstroPy
153  test_planet(venus,
154  CESkyCoord(CEAngle::Deg(183.8496760951891247),
155  CEAngle::Deg(1.8722067172672485),
157  CESkyCoord(CEAngle::Deg(242.8427330437376099),
158  CEAngle::Deg(43.8958267247570220),
160  {-0.72543765, -0.04893718, 0.02371176},
161  {0.00080326, -0.01849847, -0.00837267});
162 
163  return pass();
164 }
165 
166 
167 /**********************************************************************/
172 {
173  // Create the planet object
174  CEPlanet earth = CEPlanet::Earth();
175 
176  // Update the coordinates
178 
179  // Run the actual test with values derived from AstroPy
180  test_planet(earth,
181  CESkyCoord(CEAngle::Deg(101.7655134417398273),
182  CEAngle::Deg(23.0103434895188776),
184  CESkyCoord(CEAngle::Deg(89.9999999991386233),
185  CEAngle::Deg(179.9999111107912881),
187  {-0.18428431, 0.88477935, 0.383819},
188  {-0.01720221, -0.00290513, -0.00125952});
189 
190  return pass();
191 }
192 
193 
194 /**********************************************************************/
199 {
200  // Create the planet object
201  CEPlanet mars = CEPlanet::Mars();
202 
203  // Update the coordinates
205 
206  // Run the actual test with values derived from AstroPy
207  test_planet(mars,
208  CESkyCoord(CEAngle::Deg(359.9442433037739306),
209  CEAngle::Deg(-1.5700885004650793),
211  CESkyCoord(CEAngle::Deg(106.9869245851015762),
212  CEAngle::Deg(51.3129710961485586),
214  {1.38356924, -0.0011989, -0.0378561},
215  {0.00067763, 0.01380768, 0.00631503});
216 
217  // Test that we can return the position and velocity values
218  std::vector<double> pos = mars.PositionICRS();
219  std::vector<double> vel = mars.VelocityICRS();
220  test_double(mars.GetXICRS(), pos[0], __func__, __LINE__);
221  test_double(mars.GetYICRS(), pos[1], __func__, __LINE__);
222  test_double(mars.GetZICRS(), pos[2], __func__, __LINE__);
223  test_double(mars.GetVxICRS(), vel[0], __func__, __LINE__);
224  test_double(mars.GetVyICRS(), vel[1], __func__, __LINE__);
225  test_double(mars.GetVzICRS(), vel[2], __func__, __LINE__);
226 
227 
228  // Run the test now for the JPL algorithm. Note that the precision is going
229  // to be much worse, so we only check for rough agreement (within 0.1 deg).
230  CEPlanet mars2 = CEPlanet::Mars();
231  mars2.SetAlgorithm(CEPlanetAlgo::JPL);
232  mars2.UpdateCoordinates(base_date_.JD());
233  CEAngle angsep = mars.AngularSeparation(mars2);
234  test_lessthan(angsep.Deg(), 0.1, __func__, __LINE__);
235 
236  return pass();
237 }
238 
239 
240 /**********************************************************************/
245 {
246  // Create the planet object
247  CEPlanet jupiter = CEPlanet::Jupiter();
248 
249  // Update the coordinates
250  jupiter.UpdateCoordinates(base_date_.JD());
251 
252  // Run the actual test with values derived from AstroPy
253  test_planet(jupiter,
254  CESkyCoord(CEAngle::Deg(34.3822629405528843),
255  CEAngle::Deg(12.5158780867456674),
257  CESkyCoord(CEAngle::Deg(81.1700048918871886),
258  CEAngle::Deg(103.2491992750974532),
260  {3.99442023, 2.7334608, 1.07451899},
261  {-0.00455544, 0.00587705, 0.00263009});
262 
263  return pass();
264 }
265 
266 
267 /**********************************************************************/
272 {
273  // Create the planet object
274  CEPlanet saturn = CEPlanet::Saturn();
275 
276  // Update the coordinates
277  saturn.UpdateCoordinates(base_date_.JD());
278 
279  // Run the actual test with values derived from AstroPy
280  test_planet(saturn,
281  CESkyCoord(CEAngle::Deg(43.9734819060061355),
282  CEAngle::Deg(14.3451097907081060),
284  CESkyCoord(CEAngle::Deg(75.7370064470109696),
285  CEAngle::Deg(117.5753277202141760),
287  {6.39746262, 6.17262105, 2.2735304},
288  {-0.00429156, 0.00350834, 0.00163369});
289 
290  return pass();
291 }
292 
293 
294 /**********************************************************************/
299 {
300  // Create the planet object
301  CEPlanet uranus = CEPlanet::Uranus();
302 
303  // Update the coordinates
304  uranus.UpdateCoordinates(base_date_.JD());
305 
306  // Run the actual test with values derived from AstroPy
307  test_planet(uranus,
308  CESkyCoord(CEAngle::Deg(319.0662412483117691),
309  CEAngle::Deg(-16.5756054581048744),
311  CESkyCoord(CEAngle::Deg(116.9533213948411401),
312  CEAngle::Deg(40.2264403677069211),
314  {14.42492523, -12.50957402, -5.6830779},
315  {0.00269067, 0.00244776, 0.00103396});
316 
317  return pass();
318 }
319 
320 
321 /**********************************************************************/
326 {
327  // Create the planet object
328  CEPlanet neptune = CEPlanet::Neptune();
329 
330  // Update the coordinates
331  neptune.UpdateCoordinates(base_date_.JD());
332 
333  // Run the actual test with values derived from AstroPy
334  test_planet(neptune,
335  CESkyCoord(CEAngle::Deg(306.1731137010386306),
336  CEAngle::Deg(-19.0396553491478997),
338  CESkyCoord(CEAngle::Deg(129.5362400208609870),
339  CEAngle::Deg(31.1291412849027509),
341  {16.80489053, -22.98266171, -9.82533257},
342  {0.00258607, 0.00165554, 0.00061312});
343 
344  return pass();
345 }
346 
347 
348 /**********************************************************************/
357 bool test_CEPlanet::test_planet(const CEPlanet& test_planet,
358  const CESkyCoord& true_icrs,
359  const CESkyCoord& true_obs,
360  const std::vector<double>& true_pos,
361  const std::vector<double>& true_vel)
362 {
363  // Define a name so that we actually know which planet we're testing
364  std::string func_name = std::string(__func__) + " (" + test_planet.Name() + ")";
365 
366  // Test ICRS coordinates
367  CESkyCoord icrs_coords(test_planet.XCoord(),
368  test_planet.YCoord(),
370  if (!test(icrs_coords == true_icrs, func_name, __LINE__)) {
371  // Print information about the coordinates
372  std::printf(" %s %s", icrs_coords.print().c_str(), true_icrs.print().c_str());
373  std::printf(" X-diff: %f arcsec\n", (icrs_coords.XCoord().Deg()-true_icrs.XCoord().Deg())*3600.0);
374  std::printf(" Y-diff: %f arcsec\n", (icrs_coords.YCoord().Deg()-true_icrs.YCoord().Deg())*3600.0);
375  }
376  std::printf(" AngSep ICRS: %e arcsec\n", icrs_coords.AngularSeparation(true_icrs).Deg()*3600.0);
377 
378  // Test observed coordinates
379  CESkyCoord obs_coords = test_planet.ObservedCoords(base_date_,
381 
382  if (!test(obs_coords == true_obs, func_name, __LINE__)) {
383  // Print information about the coordinates
384  std::printf(" %s %s", obs_coords.print().c_str(), true_obs.print().c_str());
385  std::printf(" X-diff: %f arcsec\n", (obs_coords.XCoord().Deg()-true_obs.XCoord().Deg())*3600.0);
386  std::printf(" Y-diff: %f arcsec\n", (obs_coords.YCoord().Deg()-true_obs.YCoord().Deg())*3600.0);
387  }
388  std::printf(" AngSep Obs : %e arcsec\n", obs_coords.AngularSeparation(true_obs).Deg()*3600.0);
389 
390  // Update the tolerance
391  double tol_old = DblTol();
392  SetDblTol(1.0e-5);
393 
394  // Test the computed ICRS position values
395  std::vector<double> test_pos = test_planet.PositionICRS();
396  if (!test_vect(test_pos, true_pos, func_name, __LINE__)) {
397  double dist(0.0);
398  for (int i=0; i<test_pos.size(); i++) {
399  double offset = (test_pos[i]-true_pos[i]);
400  std::printf(" test: %f AU | expected: %f AU (rel. diff: %e AU)\n",
401  test_pos[i], true_pos[i], offset/true_pos[i]);
402  dist += (offset * offset);
403  }
404  std::printf(" Offset: %e AU\n", std::sqrt(dist));
405  }
406 
407  // Test the computed ICRS velocities
408  std::vector<double> test_vel = test_planet.VelocityICRS();
409  if (!test_vect(test_vel, true_vel, func_name, __LINE__)) {
410  double dist(0.0);
411  for (int i=0; i<test_vel.size(); i++) {
412  double offset = (test_vel[i]-true_vel[i]);
413  std::printf(" test: %e AU/day | expected: %e AU/day (rel. diff: %e AU/day)\n",
414  test_vel[i], true_vel[i], offset/true_vel[i]);
415  dist += (offset * offset);
416  }
417  std::printf(" Offset: %e AU/day\n", std::sqrt(dist));
418  }
419 
420  // Reset tolerance
421  SetDblTol(tol_old);
422 
423  return pass();
424 }
425 
426 
427 /**********************************************************************/
433 {
434  test_mercury();
435  test_venus();
436  test_earth();
437  test_mars();
438  test_jupiter();
439  test_saturn();
440  test_uranus();
441  test_neptune();
442 
443  return pass();
444 }
445 
446 
447 /**********************************************************************/
450 int main(int argc, char** argv)
451 {
452  test_CEPlanet tester;
453  return (!tester.runtests());
454 }
test_CEPlanet::base_observer_
CEObserver base_observer_
Definition: test_CEPlanet.h:79
test_CEPlanet::test_mars
virtual bool test_mars(void)
Test Mars.
Definition: test_CEPlanet.cpp:198
CESkyCoordType::ICRS
RA, Dec (referenced at the barycenter of the solarsystem)
CEPlanet::Uranus
static CEPlanet Uranus()
Returns an object representing Uranus.
Definition: CEPlanet.cpp:391
test_CEPlanet::base_date_
CEDate base_date_
Definition: test_CEPlanet.h:80
CENamespace.h
main
int main(int argc, char **argv)
Main method that actually runs the tests.
Definition: test_CEPlanet.cpp:450
CEDate
Definition: CEDate.h:43
CEPlanet::Mercury
static CEPlanet Mercury()
Returns an object representing Mercury.
Definition: CEPlanet.cpp:211
test_CEPlanet::test_jupiter
virtual bool test_jupiter(void)
Test Jupiter.
Definition: test_CEPlanet.cpp:244
test_CEPlanet::test_earth
virtual bool test_earth(void)
Test Earth.
Definition: test_CEPlanet.cpp:171
CEPlanet::Mars
static CEPlanet Mars()
Returns an object representing Mars.
Definition: CEPlanet.cpp:308
CEPlanet::Neptune
static CEPlanet Neptune()
Returns an object representing Neptune.
Definition: CEPlanet.cpp:419
CESkyCoord
Definition: CESkyCoord.h:49
test_CEPlanet::test_uranus
virtual bool test_uranus(void)
Test Uranus.
Definition: test_CEPlanet.cpp:298
test_CEPlanet::test_CEPlanet
test_CEPlanet()
Default constructor.
Definition: test_CEPlanet.cpp:29
test_CEPlanet::~test_CEPlanet
virtual ~test_CEPlanet()
Destructor.
Definition: test_CEPlanet.cpp:52
test_CEPlanet::test_construct
virtual bool test_construct(void)
Test constructors.
Definition: test_CEPlanet.cpp:78
CESkyCoord::AngularSeparation
virtual CEAngle AngularSeparation(const CESkyCoord &coords) const
Get the angular separation between the coordinates represented by this object and another coordinate ...
Definition: CESkyCoord.cpp:119
CEBody::GetCoordinates
virtual CESkyCoord GetCoordinates(const CEDate &date=CEDate::CurrentJD()) const
Return the ICRS coordinates associated with this object.
Definition: CEBody.h:121
test_CEPlanet::test_saturn
virtual bool test_saturn(void)
Test Saturn.
Definition: test_CEPlanet.cpp:271
CEObserver::SetRelativeHumidity
void SetRelativeHumidity(const double &humidity=0.0)
Set the observer's relative humidity.
Definition: CEObserver.h:349
CEPlanet::GetXICRS
double GetXICRS()
Return X distance from solar system barycenter (AU)
Definition: CEPlanet.h:320
CEObserver::SetWavelength_um
void SetWavelength_um(const double &new_wavelength_um)
Set the observer's observing wavelength (micrometers)
Definition: CEObserver.h:393
CEPlanet::Earth
static CEPlanet Earth()
Returns an object representing Earth.
Definition: CEPlanet.cpp:261
CEObserver::SetPressure_hPa
void SetPressure_hPa(const double &pressure=CppEphem::EstimatePressure_hPa(CppEphem::SeaLevelTemp_C()))
Set the observer's pressure.
Definition: CEObserver.h:338
CEPlanet::VelocityICRS
std::vector< double > VelocityICRS(void) const
Vector of velocities relative to solar system barycenter (AU/day)
Definition: CEPlanet.h:390
CEAngle
Definition: CEAngle.h:38
test_CEPlanet.h
JPL
Use Keplerian algorithm outlined by JPL.
Definition: CEPlanet.h:49
CEObserver::SetTemperature_C
void SetTemperature_C(const double &temp_C=CppEphem::SeaLevelTemp_C())
Set the observer's temperature (Celsius)
Definition: CEObserver.h:360
CppEphem::julian_date_J2000
double julian_date_J2000()
Julian Date corresponding to J2000.
Definition: CENamespace.h:67
CESkyCoord::GetCoordSystem
CESkyCoordType GetCoordSystem(void) const
Return coordinate system.
Definition: CESkyCoord.h:251
CEPlanet::UpdateCoordinates
virtual void UpdateCoordinates(double new_date=-1.0e30) const
Recomputes the coordinates of the planet if the date has changed from the last time the coordinates w...
Definition: CEPlanet.cpp:677
CEAngleType::DEGREES
JD
Julian Date.
Definition: CEDate.h:56
CEBody::Name
std::string Name(void) const
Get the name of this object.
Definition: CEBody.h:95
CEPlanet::Venus
static CEPlanet Venus()
Returns an object representing Venus.
Definition: CEPlanet.cpp:236
CEPlanet::YCoordinate_Deg
virtual double YCoordinate_Deg(double new_date=-1.0e30) const
Definition: CEPlanet.h:309
test_CEPlanet::runtests
virtual bool runtests()
Run tests.
Definition: test_CEPlanet.cpp:61
test_CEPlanet::test_neptune
virtual bool test_neptune(void)
Test Neptune.
Definition: test_CEPlanet.cpp:325
CEPlanet
Definition: CEPlanet.h:35
CEPlanet::GetVxICRS
double GetVxICRS()
X velocity relative to solar system barycenter (AU/day)
Definition: CEPlanet.h:360
test_CEPlanet::test_mercury
virtual bool test_mercury(void)
Test Mercury.
Definition: test_CEPlanet.cpp:117
CEPlanet::GetYICRS
double GetYICRS()
Y distance from solar system barycenter (AU)
Definition: CEPlanet.h:330
CESkyCoord::print
std::string print(void) const
Generate a message string that specifies the information about this coordinate.
Definition: CESkyCoord.cpp:1026
CEPlanet::Saturn
static CEPlanet Saturn()
Returns an object representing Saturn.
Definition: CEPlanet.cpp:363
test_CEPlanet
Definition: test_CEPlanet.h:29
CEObserver
Definition: CEObserver.h:28
CEAngle::Deg
double Deg(void) const
Return angle in degrees as a double.
Definition: CEAngle.cpp:328
test_CEPlanet::test_venus
virtual bool test_venus(void)
Test Venus.
Definition: test_CEPlanet.cpp:144
CppEphem::CorrectionsInterp
void CorrectionsInterp(bool set_interp)
Set the corrections object to use interpolation.
Definition: CENamespace.cpp:99
test_CEPlanet::test_Planets
bool test_Planets(void)
Compare.
Definition: test_CEPlanet.cpp:432
CEDate::SetReturnType
void SetReturnType(CEDateType return_type)
Set the return type from the overloaded 'operator double'.
Definition: CEDate.h:247
CEPlanet::PositionICRS
std::vector< double > PositionICRS(void) const
Z distance from solar system barycenter (AU)
Definition: CEPlanet.h:350
CEDate::JD
virtual double JD() const
Get the Julian date represented by this object.
Definition: CEDate.h:150
CEPlanet::GetVyICRS
double GetVyICRS()
Y velocity relative to solar system barycenter (AU/day)
Definition: CEPlanet.h:370
CEPlanet::Jupiter
static CEPlanet Jupiter()
Returns an object representing Jupiter.
Definition: CEPlanet.cpp:335
CEPlanet::XCoordinate_Deg
virtual double XCoordinate_Deg(double new_date=-1.0e30) const
Definition: CEPlanet.h:289
CEAngle::Deg
static CEAngle Deg(const double &angle)
Return angle (radians) constructed from a degree angle.
Definition: CEAngle.cpp:317
CEPlanet::GetVzICRS
double GetVzICRS()
Z velocity relative to solar system barycenter (AU/day)
Definition: CEPlanet.h:380
CESkyCoord::XCoord
virtual CEAngle XCoord(const CEDate &jd=CppEphem::julian_date_J2000()) const
Return x coordinate at given Julian date.
Definition: CESkyCoord.h:227
test_CEPlanet::test_planet
virtual bool test_planet(const CEPlanet &test_planet, const CESkyCoord &true_icrs, const CESkyCoord &true_obs, const std::vector< double > &true_pos, const std::vector< double > &true_vel)
Test parameters for a given planet.
Definition: test_CEPlanet.cpp:357
CESkyCoordType::OBSERVED
Azimuth, Zenith (requires additional observer information)
CESkyCoord::YCoord
virtual CEAngle YCoord(const CEDate &jd=CppEphem::julian_date_J2000()) const
Return y coordinate at given Julian date.
Definition: CESkyCoord.h:240
CEPlanet::GetZICRS
double GetZICRS()
Z distance from solar system barycenter (AU)
Definition: CEPlanet.h:340