Planetary Orbital Evolution due to Tides
Orbital evolution of two objects experiencing tides
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