Planetary Orbital Evolution due to Tides
Orbital evolution of two objects experiencing tides
•All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Modules Pages
testDifferentialEquations.cpp
Go to the documentation of this file.
1 
8 #include "testDifferentialEquations.h"
9 
10 namespace Evolve {
11 
12  void test_DifferentialEquations::output_rates(
13  const alglib::real_1d_array &eccentricities,
14  const alglib::real_1d_array &expected_semimajor_rate,
15  const alglib::real_1d_array &predicted_semimajor_rate,
16  const alglib::real_1d_array &expected_eccentricity_rate,
17  const alglib::real_1d_array &predicted_eccentricity_rate
18  ) const
19  {
20  std::cout << std::setw(25) << "eccentricity"
21  << std::setw(25) << "semimajor_rate"
22  << std::setw(25) << "expected_semimajor_rate"
23  << std::setw(25) << "eccentricity_rate"
24  << std::setw(25) << "expected_ecc_rate"
25  << std::endl;
26  for(
27  unsigned e_index = 0;
28  e_index < eccentricities.length();
29  ++e_index
30  ) {
31  std::cout << std::setw(25) << eccentricities[e_index]
32  << std::setw(25)
33  << predicted_semimajor_rate[e_index]
34  << std::setw(25)
35  << expected_semimajor_rate[e_index]
36  << std::setw(25)
37  << predicted_eccentricity_rate[e_index]
38  << std::setw(25)
39  << expected_eccentricity_rate[e_index]
40  << std::endl;
41  }
42 
43 
44  }
45 
46  double test_DifferentialEquations::zahn77_semimajor_rate_coef(
47  int orbital_frequency_multiplier,
48  int spin_frequency_multiplier,
49  double eccentricity
50  ) const
51  {
52 
53  double e2 = std::pow(eccentricity, 2);
54 
55  if(orbital_frequency_multiplier == 2 && spin_frequency_multiplier == 2)
56  return 1.0 - 5.0 * e2;
57 
58  if(orbital_frequency_multiplier == 1 && spin_frequency_multiplier == 0)
59  return 3.0 / 4.0 * e2;
60 
61  if(orbital_frequency_multiplier == 1 && spin_frequency_multiplier == 2)
62  return 1.0 / 8.0 * e2;
63 
64  if(orbital_frequency_multiplier == 3 && spin_frequency_multiplier == 2)
65  return 147.0 / 8.0 * e2;
66 
67  return 0.0;
68  }
69 
70  double test_DifferentialEquations::zahn77_eccentricity_rate_coef(
71  int orbital_frequency_multiplier,
72  int spin_frequency_multiplier
73  ) const
74  {
75 
76  if(orbital_frequency_multiplier == 2 && spin_frequency_multiplier == 2)
77  return -1.0;
78 
79  if(orbital_frequency_multiplier == 1 && spin_frequency_multiplier == 0)
80  return 1.5;
81 
82  if(orbital_frequency_multiplier == 1 && spin_frequency_multiplier == 2)
83  return -1.0 / 4.0;
84 
85  if(orbital_frequency_multiplier == 3 && spin_frequency_multiplier == 2)
86  return 49.0 / 4.0;
87 
88  return 0.0;
89  }
90 
91  void test_DifferentialEquations::check_agreement(
92  const alglib::real_1d_array& x,
93  const alglib::real_1d_array& y1,
94  const alglib::real_1d_array& y2,
95  unsigned agreement_order,
96  unsigned max_order,
97  const std::string &description
98  )
99  {
100  alglib::real_1d_array difference;
101  difference.setlength(x.length());
102  double min_difference = Core::Inf,
103  max_difference = -Core::Inf,
104  y_scale = 0;
105  for(unsigned i = 0; i < x.length(); ++i) {
106  difference[i] = y2[i] - y1[i];
107  min_difference = std::min(difference[i], min_difference);
108  max_difference = std::max(difference[i], max_difference);
109  y_scale = std::max(std::abs(y1[i]), std::abs(y2[i]));
110  }
111 
112  if(
113  std::max(std::abs(max_difference), std::abs(min_difference))
114  <
115  1e-13 * y_scale
116  )
117  return;
118  alglib::ae_int_t fit_info;
119  alglib::barycentricinterpolant interp;
120  alglib::polynomialfitreport fit_report;
121  alglib::polynomialfit(x,
122  difference,
123  max_order + 1,
124  fit_info,
125  interp,
126  fit_report);
127 
128  double max_allowed_residual = 1e-12 * (max_difference - min_difference);
129 
130  std::ostringstream message;
131  message << description
132  << " too large residual from fit to difference: "
133  << fit_report.maxerror
134  << " > "
135  << max_allowed_residual;
136 
137  TEST_ASSERT_MSG(fit_report.maxerror <= max_allowed_residual,
138  message.str().c_str());
139 
140  alglib::real_1d_array coefficients;
141 
142  alglib::polynomialbar2pow(interp, coefficients);
143 
144  double max_allowed_term = 0.0,
145  max_abs_x = std::max(std::abs(x[0]), x[x.length() - 1]);
146  for(unsigned power = 0; power < coefficients.length(); ++power)
147  max_allowed_term = std::max(
148  max_allowed_term,
149  std::abs(coefficients[power]) * std::pow(max_abs_x, power)
150  );
151  max_allowed_term *= 1e-8;
152 
153  for(unsigned power = 0; power <= agreement_order; ++power) {
154  double max_abs_term_value = (std::abs(coefficients[power])
155  *
156  std::pow(max_abs_x, power));
157  message.str("");
158  message << description
159  << " x^" << power << " coefficient too large: |"
160  << coefficients[power]
161  << " * "
162  << max_abs_x << "^" << power
163  << "| = "
164  << max_abs_term_value
165  << " > "
166  << max_allowed_term;
167 
168  TEST_ASSERT_MSG(max_abs_term_value <= max_allowed_term,
169  message.str().c_str());
170  }
171  }
172 
173  test_DifferentialEquations::test_DifferentialEquations()
174  {
175  TEST_ADD(test_DifferentialEquations::test_aligned_equations);
176  TEST_ADD(test_DifferentialEquations::test_error_estimate);
177  }
178 
179  void test_DifferentialEquations::test_aligned_equations()
180  {
181  const double MAX_ECCENTRICITY = 0.5;
182  const unsigned NUM_ECCENTRICITIES = 100;
183  const double ECCENTRICITY_STEP = (
184  MAX_ECCENTRICITY / (NUM_ECCENTRICITIES - 1)
185  );
186 
187  alglib::real_1d_array eccentricities,
188  expected_semimajor_rate,
189  predicted_semimajor_rate,
190  expected_eccentricity_rate,
191  predicted_eccentricity_rate;
192  eccentricities.setlength(NUM_ECCENTRICITIES);
193  expected_semimajor_rate.setlength(NUM_ECCENTRICITIES);
194  predicted_semimajor_rate.setlength(NUM_ECCENTRICITIES);
195  expected_eccentricity_rate.setlength(NUM_ECCENTRICITIES);
196  predicted_eccentricity_rate.setlength(NUM_ECCENTRICITIES);
197 
198  eccentricities[0] = 0.0;
199  for(unsigned e_index = 1; e_index < NUM_ECCENTRICITIES; ++e_index)
200  eccentricities[e_index] = (eccentricities[e_index - 1]
201  +
202  ECCENTRICITY_STEP);
203 
204  std::valarray<double> one_poly(1.0, 1);
205  for(
206  int orbital_frequency_multiplier = 0;
207  orbital_frequency_multiplier <= 4;
208  ++orbital_frequency_multiplier
209  ) {
210  for(
211  int spin_frequency_multiplier = 0;
212  spin_frequency_multiplier <= 4;
213  ++spin_frequency_multiplier
214  ) {
215  SingleTidalTermBody primary(0.0,
216  1.0,
217  orbital_frequency_multiplier,
218  spin_frequency_multiplier,
219  1.0,
220  one_poly,
221  one_poly,
222  one_poly);
223  Planet::Planet secondary(1.0, 1.0);
224 
225  BinarySystem binary(primary, secondary);
226 
227  binary.change_e_order(2);
228 
229  std::valarray<double> parameters(0.0, 7);
230  std::valarray<double> evolution_rates(parameters.size());
231  parameters[0] = 10.0;
232  parameters[1] = 0.0;
233  binary.configure(true,
234  0.0,
235  &(parameters[0]),
236  Core::BINARY);
237 
238  double common_rate_factor = (
239  -2.4 * M_PI
240  *
241  (primary.mass() + secondary.mass()) / primary.mass()
242  *
243  Core::AstroConst::G
244  *
245  Core::AstroConst::solar_mass * secondary.mass()
246  /
247  std::pow(Core::AstroConst::solar_radius, 3)
248  ) * Core::AstroConst::Gyr;
249 
250  for(double semimajor = 2.0; semimajor <= 10.0; semimajor += 1.0) {
251 
252  double mean_motion = std::sqrt(
253  Core::AstroConst::G * Core::AstroConst::solar_mass * (
254  primary.mass()
255  +
256  secondary.mass()
257  )
258  /
259  std::pow(semimajor * Core::AstroConst::solar_radius, 3)
260  );
261  double rate_factor = (
262  common_rate_factor
263  /
264  (mean_motion * std::pow(semimajor, 8))
265  );
266  for(
267  double primary_spin_frequency = 0.0;
268  primary_spin_frequency <= 5.0 * mean_motion;
269  primary_spin_frequency += 0.2 * mean_motion
270  ) {
271 
272  for(
273  unsigned e_index = 0;
274  e_index < NUM_ECCENTRICITIES;
275  ++e_index
276  )
277  {
278  double age = 1.0,
279  test_e = eccentricities[e_index];
280  parameters[0] = std::pow(semimajor, 6.5);
281  parameters[1] = test_e;
282  binary.differential_equations(age,
283  &(parameters[0]),
284  Core::BINARY,
285  &(evolution_rates[0]));
286  expected_semimajor_rate[e_index] = (
287  6.5 * std::pow(semimajor, 6.5)
288  *
289  rate_factor
290  *
291  zahn77_semimajor_rate_coef(
292  orbital_frequency_multiplier,
293  spin_frequency_multiplier,
294  test_e
295  )
296  );
297  expected_eccentricity_rate[e_index] = (
298  test_e
299  *
300  rate_factor / 4.0
301  *
302  zahn77_eccentricity_rate_coef(
303  orbital_frequency_multiplier,
304  spin_frequency_multiplier
305  )
306  );
307  predicted_semimajor_rate[e_index] = evolution_rates[0];
308  predicted_eccentricity_rate[e_index] =
309  evolution_rates[1];
310  }
311 /* output_rates(eccentricities,
312  expected_semimajor_rate,
313  predicted_semimajor_rate,
314  expected_eccentricity_rate,
315  predicted_eccentricity_rate);*/
316  std::ostringstream message;
317  message << "orbital freq. mult. = "
318  << orbital_frequency_multiplier
319  << ", spin freq. mult. = "
320  << spin_frequency_multiplier
321  << " for a = " << semimajor
322  << ", W* = " << primary_spin_frequency;
323  check_agreement(
324  eccentricities,
325  expected_semimajor_rate,
326  predicted_semimajor_rate,
327  2,
328  4,
329  (
330  "Checking semimajor rate against Zahn (1977)"
331  +
332  message.str()
333  )
334  );
335  check_agreement(
336  eccentricities,
337  expected_eccentricity_rate,
338  predicted_eccentricity_rate,
339  2,
340  20,
341  (
342  "Checking eccentricity rate against Zahn (1977)"
343  +
344  message.str()
345  )
346  );
347  }
348  }
349  }
350  }
351  }
352 
353  void test_DifferentialEquations::test_error_estimate()
354  {
355  const double MAX_ECCENTRICITY = 0.9;
356  const unsigned NUM_ECCENTRICITIES = 100;
357  const double ECCENTRICITY_STEP = (
358  MAX_ECCENTRICITY / (NUM_ECCENTRICITIES - 1)
359  );
360  const double a = 10.0;
361  const double age = 1.0;
362 
363  StellarEvolution::MockStellarEvolution *no_evol =
364  StellarEvolution::make_no_evolution();
365 
366  Star::InterpolatedEvolutionStar *star = make_const_lag_star(
367  *no_evol,
368  1.0,
369  1.0,
370  1.0,
371  1.0
372  );
373 
374  Planet::Planet planet(1.0, 1.0);
375 
376  BinarySystem binary(*star, planet);
377 
378  std::valarray<double> parameters(0.0, 10),
379  rough_rates(parameters.size()),
380  fine_rates(parameters.size()),
381  rate_errors(parameters.size());
382 
383  parameters[0] = a;
384 
385  binary.configure(true,
386  (MIN_AGE + MAX_AGE) / 2.0,
387  &(parameters[0]),
388  Core::BINARY);
389  star->detect_saturation();
390 
391  std::cout << std::setw(25) << "e_order"
392  << std::setw(25) << "e"
393  << std::setw(25) << "rough_a_rate"
394  << std::setw(25) << "fine_a_rate"
395  << std::setw(25) << "a_rate_error"
396  << std::setw(25) << "rough_e_rate"
397  << std::setw(25) << "fine_e_rate"
398  << std::setw(25) << "e_rate_error"
399  << std::setw(25) << "rough_inc_rate"
400  << std::setw(25) << "fine_inc_rate"
401  << std::setw(25) << "inc_rate_error"
402  << std::setw(25) << "rough_spin_rate"
403  << std::setw(25) << "fine_spin_rate"
404  << std::setw(25) << "spin_rate_error"
405 
406  << std::endl;
407 
408 
409  for(unsigned e_order = 5; e_order <= 100; e_order+=5) {
410  for(double e = 0.0; e <= MAX_ECCENTRICITY; e += ECCENTRICITY_STEP) {
411  parameters[1] = e;
412 
413  star->zone(0).change_e_order(e_order, binary, true, 0);
414  binary.differential_equations(age,
415  &(parameters[0]),
416  Core::BINARY,
417  &(rough_rates[0]));
418  binary.differential_equations(age,
419  &(parameters[0]),
420  Core::BINARY,
421  &(rate_errors[0]),
422  true);
423 
424  star->zone(0).change_e_order(2 * e_order, binary, true, 0);
425  binary.differential_equations(age,
426  &(parameters[0]),
427  Core::BINARY,
428  &(fine_rates[0]));
429 
430  std::cout << std::setw(25) << e_order
431  << std::setw(25) << e
432  << std::setw(25) << rough_rates[0]
433  << std::setw(25) << fine_rates[0]
434  << std::setw(25) << rate_errors[0]
435  << std::setw(25) << rough_rates[1]
436  << std::setw(25) << fine_rates[1]
437  << std::setw(25) << rate_errors[1]
438  << std::setw(25) << rough_rates[2]
439  << std::setw(25) << fine_rates[2]
440  << std::setw(25) << rate_errors[2]
441  << std::setw(25) << rough_rates[7]
442  << std::setw(25) << fine_rates[7]
443  << std::setw(25) << rate_errors[7]
444  << std::endl;
445  }
446  }
447 
448  delete star;
449  delete no_evol;
450  }
451 
452 } //End Envolve namespace.
StellarEvolution::MockStellarEvolution
Implements a StellarEvolution based on polynomial evolution quantities.
Definition: PolynomialEvolution.h:126
MIN_AGE
const double MIN_AGE
Most tests start at this age in Gyr.
Definition: Common.h:54
Evolve::BinarySystem::differential_equations
int differential_equations(double age, const double *parameters, Core::EvolModeType evolution_mode, double *differential_equations, bool expansion_error=false)
The differential equation and jacobian for the evolution of the system.
Definition: BinarySystem.cpp:2106
Star::SaturatingSkumanichWindBody::detect_saturation
void detect_saturation()
Sets the saturation based on the currently configured spin frequency.
Definition: SaturatingSkumanichWindBody.h:65
Evolve::test_DifferentialEquations::test_DifferentialEquations
test_DifferentialEquations()
Create the test suite.
Definition: testDifferentialEquations.cpp:173
Evolve::SingleTidalTermBody
A skumanich wind body with a single zone dissipative to only a single tidal term. ...
Definition: SingleTidalTermBody.h:21
Evolve::test_DifferentialEquations::zahn77_semimajor_rate_coef
double zahn77_semimajor_rate_coef(int orbital_frequency_multiplier, int spin_frequency_multiplier, double eccentricity) const
The Zahn (1977) coefficient for the semimajor evolution rate that corresponds to the given tidal term...
Definition: testDifferentialEquations.cpp:46
Evolve::test_DifferentialEquations::test_aligned_equations
void test_aligned_equations()
Test the a & e differential equations for aligned orbit.
Definition: testDifferentialEquations.cpp:179
Planet::Planet
Single zone non-evolving planets with huge dissipation, so they always remain locked to the disk...
Definition: Planet.h:21
Star::InterpolatedEvolutionStar::zone
Evolve::DissipatingZone & zone(unsigned zone_index)
See Evolve::DissipatingBody::zone().
Definition: EvolvingStar.h:85
MAX_AGE
const double MAX_AGE
Most tests end at this age in Gyr.
Definition: Common.h:57
Evolve
Orientations of zones of bodies in a binary system.
Definition: ConstSolutionIterator.cpp:3
Evolve::DissipatingBody::mass
double mass() const
The mass of the body (constant with age).
Definition: DissipatingBody.h:534
Core::AstroConst::solar_radius
const double solar_radius
Radius of the sun [m].
Definition: AstronomicalConstants.h:22
Evolve::test_DifferentialEquations::check_agreement
void check_agreement(const alglib::real_1d_array &x, const alglib::real_1d_array &y1, const alglib::real_1d_array &y2, unsigned agreement_order, unsigned max_order, const std::string &description)
Check if two given dependencies on x agree up to a given order in x.
Definition: testDifferentialEquations.cpp:91
Core::AstroConst::Gyr
const double Gyr
Gyr [s].
Definition: AstronomicalConstants.h:24
Core::AstroConst::solar_mass
const double solar_mass
Mass of the sun [kg].
Definition: AstronomicalConstants.h:21
Evolve::test_DifferentialEquations::zahn77_eccentricity_rate_coef
double zahn77_eccentricity_rate_coef(int orbital_frequency_multiplier, int spin_frequency_multiplier) const
Same as zahn77_semimajor_rate_coef() but for the eccentricity evolution.
Definition: testDifferentialEquations.cpp:70
testDifferentialEquations.h
Unit tests that check the differential equations for eccentricity and semimajor against analytic expr...
Evolve::test_DifferentialEquations::test_error_estimate
void test_error_estimate()
Test the error estimate of the differential equations.
Definition: testDifferentialEquations.cpp:353
Evolve::BinarySystem::change_e_order
virtual void change_e_order(unsigned new_e_order)
Change the eccentricity expansion order for all dissipative zones.
Definition: BinarySystem.h:1060
Evolve::DissipatingZone::change_e_order
virtual void change_e_order(unsigned new_e_order, BinarySystem &system, bool primary, unsigned zone_index)
Changes the order of the eccentricity expansion performed.
Definition: DissipatingZone.cpp:902
make_const_lag_star
Star::InterpolatedEvolutionStar * make_const_lag_star(const StellarEvolution::Interpolator &evolution, double wind_strength, double wind_sat_freq, double coupling_timescale, double phase_lag)
Create a star with the given parameters with a constant phase lag.
Definition: MakeStar.cpp:10
Evolve::BinarySystem::configure
virtual int configure(bool initialize, double age, double semimajor, double eccentricity, const double *spin_angmom, const double *inclination, const double *periapsis, Core::EvolModeType evolution_mode)
Sets the current state of the system.
Definition: BinarySystem.cpp:1912
Evolve::BinarySystem
Describes a system of two bodies orbiting each other.
Definition: BinarySystem.h:56
Evolve::test_DifferentialEquations::output_rates
void output_rates(const alglib::real_1d_array &eccentricities, const alglib::real_1d_array &expected_semimajor_rate, const alglib::real_1d_array &predicted_semimajor_rate, const alglib::real_1d_array &expected_eccentricity_rate, const alglib::real_1d_array &predicted_eccentricity_rate) const
Output the predicted/expected semiamjor and eccentricity evolution rates to stdout.
Definition: testDifferentialEquations.cpp:12
Core::AstroConst::G
const double G
Gravitational constant in SI.
Definition: AstronomicalConstants.h:26