Planetary Orbital Evolution due to Tides
Orbital evolution of two objects experiencing tides
testOrbitSolver.cpp
Go to the documentation of this file.
1 
9 #include "testOrbitSolver.h"
10 
11 namespace Evolve {
12 
13  const double TSTART = 2.0 * MIN_AGE;
14  const double zero = 0.0;
15  const double one = 1.0;
16 
17  const double Mjup_to_Msun = (Core::AstroConst::jupiter_mass
18  /
19  Core::AstroConst::solar_mass);
20 
21  const double Rjup_to_Rsun = (Core::AstroConst::jupiter_radius
22  /
23  Core::AstroConst::solar_radius);
24 
25  std::ostream &operator<<(std::ostream &os,
26  RealEvolutionQuantity q)
27  {
28  switch(q) {
29  case SEMIMAJOR :
30  os << "SEMIMAJOR"; break;
31  case ECCENTRICITY:
32  os << "ECCENTRICITY"; break;
33  case CONV_INCLINATION:
34  os << "CONV_INCLINATION"; break;
35  case RAD_INCLINATION:
36  os << "RAD_INCLINATION"; break;
37  case CONV_PERIAPSIS:
38  os << "CONV_PERIAPSIS"; break;
39  case RAD_PERIAPSIS:
40  os << "RAD_PERIAPSIS"; break;
41  case CONV_ANGMOM:
42  os << "CONV_ANGMOM"; break;
43  case RAD_ANGMOM:
44  os << "RAD_ANGMOM"; break;
45  case AGE:
46  os << "AGE"; break;
47  default :
48  assert(false);
49  }
50  return os;
51  }
52 
53 
54  StellarEvolution::PolynomialEvolutionQuantity nan_func(
55  std::valarray<double>(Core::NaN, 2),
56  TSTART,
57  MAX_AGE
58  );
59 
60  StellarEvolution::PolynomialEvolutionQuantity zero_func(
61  std::valarray<double>(),
62  TSTART,
63  MAX_AGE
64  );
65 
66  StellarEvolution::PolynomialEvolutionQuantity one_func(
67  std::valarray<double>(1.0, 1),
68  TSTART,
69  MAX_AGE
70  );
71 
72  StellarEvolution::PolynomialEvolutionQuantity two_func(
73  std::valarray<double>(2.0, 1),
74  TSTART,
75  MAX_AGE
76  );
77 
78  StellarEvolution::PolynomialEvolutionQuantity two_hundred_func(
79  std::valarray<double>(200.0, 1),
80  TSTART,
81  MAX_AGE
82  );
83 
84  void test_OrbitSolver::set_single_component_dissipation(
85  double min_frequency,
86  double max_frequency,
87  double decay_scale,
88  double phase_lag
89  )
90  {
91  //phase lag(min_frequency - decay_scale)
92  //=
93  //phase lag(max_frequency + decay_scale)
94  //=
95  //suppression_factor * phase_lag
96  const double suppression_factor = 0.01;
97 
98  std::vector<double> breaks(2);
99  breaks[0] = min_frequency;
100  breaks[1] = max_frequency;
101 
102  std::vector<double> powerlaw_indices(3);
103  powerlaw_indices[0] = (
104  std::log(suppression_factor)
105  /
106  std::log(1.0 - decay_scale / min_frequency)
107  );
108  powerlaw_indices[1] = 0.0;
109  powerlaw_indices[2] = (
110  std::log(suppression_factor)
111  /
112  std::log(1.0 + decay_scale / max_frequency)
113  );
114 
115  BrokenPowerlawPhaseLagZone *zone;
116 
117  if(__star) {
118  assert(__primary_planet == NULL);
119  zone = &(__star->envelope());
120  } else {
121  assert(__primary_planet);
122  zone = &(__primary_planet->zone());
123  }
124 
125  zone->setup(breaks,
126  std::vector<double>(),
127  powerlaw_indices,
128  std::vector<double>(1, 0.0),
129  phase_lag);
130  }
131 
132  void test_OrbitSolver::make_single_component_star(
133  const StellarEvolution::Interpolator &evolution,
134  double wind_strength,
135  double wind_sat_freq,
136  double coupling_timescale,
137  double min_frequency,
138  double max_frequency,
139  double decay_scale,
140  double phase_lag
141  )
142  {
143  __star = new Star::InterpolatedEvolutionStar(1.0,//mass
144  0.0,//feh
145  wind_strength,
146  wind_sat_freq,
147  coupling_timescale,
148  evolution);
149  set_single_component_dissipation(min_frequency,
150  max_frequency,
151  decay_scale,
152  phase_lag);
153 
154  }
155 
156  void test_OrbitSolver::evolve(double wdisk,
157  double tdisk,
158  double initial_a,
159  const double *initial_Lstar,
160  double initial_incl,
161  double secondary_mass,
162  double tsecondary,
163  double max_age,
164  double secondary_radius,
165  double precision,
166  double max_step_factor)
167  {
168  Evolve::DissipatingBody *primary;
169  if(__star)
170  primary = __star;
171  else
172  primary = __primary_planet;
173 
174  secondary_mass *= Mjup_to_Msun;
175  secondary_radius *= Rjup_to_Rsun;
176 
177  if(std::isnan(tsecondary)) tsecondary = tdisk;
178 
179  Planet::Planet secondary(secondary_mass,
180  (secondary_mass ? secondary_radius : 0.0));
181  secondary.configure(true, //init
182  tsecondary, //age
183  primary->mass(), //mass
184  initial_a, //semimajor
185  0.0, //eccentricity
186  &zero, //spin angmom
187  NULL, //inclination
188  NULL, //periapsis
189  false, //locked surface
190  true, //zero outer inclination
191  true);//zero outer periapsis
192 
193  __system = new Evolve::DiskBinarySystem(
194  *primary,
195  secondary,
196  initial_a, //semimajor
197  0.0, //eccentricity
198  initial_incl, //inclination
199  wdisk, //Wdisk
200  tdisk, //disk dissipation age
201  tsecondary //secondary formation age
202  );
203  double zeros[] = {0.0, 0.0};
204  if(__star) {
205  if(tdisk <= TSTART) {
206  if(tsecondary <= TSTART) {
207  __system->configure(true, //init
208  TSTART,
209  initial_a, //semimajor
210  0.0, //eccentricity
211  initial_Lstar, //spin angmom
212  zeros, //inclination
213  zeros, //periapsis
214  Core::BINARY);
215  } else {
216  __system->configure(true, //init
217  TSTART,
218  Core::NaN, //semimajor
219  Core::NaN, //eccentricity
220  initial_Lstar, //spin angmom
221  zeros, //inclination
222  zeros, //periapsis
223  Core::SINGLE);
224  }
225  } else {
226  __system->configure(true, //init
227  TSTART,
228  Core::NaN, //semimajor
229  Core::NaN, //eccentricity
230  initial_Lstar, //spin angmom
231  NULL, //inclination
232  NULL, //periapsis
233  Core::LOCKED_SURFACE_SPIN);
234  }
235  __star->detect_saturation();
236  } else {
237  double initial_inclinations[] = {initial_incl, 0.0, 0.0};
238  __system->configure(true, //init
239  tsecondary,
240  initial_a, //semimajor
241  0.0, //eccentricity
242  initial_Lstar, //spin angmom
243  initial_inclinations, //inclination
244  zeros, //periapsis
245  Core::BINARY);
246  }
247  __solver = new Evolve::OrbitSolver(max_age, precision);
248  (*__solver)(*__system,
249  (max_age - __system->age()) * max_step_factor,//time step
250  std::list<double>()); //no required ages
251  }
252 
253  std::vector< const std::list<double> *> test_OrbitSolver::get_evolution()
254  const
255  {
256  std::vector< const std::list<double> *>
257  tabulated_real_quantities(NUM_REAL_QUANTITIES);
258 
259  tabulated_real_quantities[AGE] = &(__solver->evolution_ages());
260 
261  tabulated_real_quantities[SEMIMAJOR] =
262  &(__system->semimajor_evolution());
263 
264  tabulated_real_quantities[ECCENTRICITY] =
265  &(__system->eccentricity_evolution());
266 
267  Evolve::DissipatingZone *primary_envelope;
268 
269  if(__star) {
270  tabulated_real_quantities[RAD_INCLINATION] =
271  &(__star->core().get_evolution_real(Evolve::INCLINATION));
272 
273  tabulated_real_quantities[RAD_PERIAPSIS] =
274  &(__star->core().get_evolution_real(Evolve::PERIAPSIS));
275 
276  tabulated_real_quantities[RAD_ANGMOM] = &(
277  __star->core().get_evolution_real(Evolve::ANGULAR_MOMENTUM)
278  );
279  primary_envelope = &(__star->envelope());
280  } else {
281  tabulated_real_quantities[RAD_INCLINATION] = NULL;
282  tabulated_real_quantities[RAD_PERIAPSIS] = NULL;
283  tabulated_real_quantities[RAD_ANGMOM] = NULL;
284  primary_envelope = &(__primary_planet->zone());
285  }
286 
287  tabulated_real_quantities[CONV_INCLINATION] = &(
288  primary_envelope->get_evolution_real(Evolve::INCLINATION)
289  );
290 
291  tabulated_real_quantities[CONV_PERIAPSIS] = &(
292  primary_envelope->get_evolution_real(Evolve::PERIAPSIS)
293  );
294 
295  tabulated_real_quantities[CONV_ANGMOM] = &(
296  primary_envelope->get_evolution_real(
297  Evolve::ANGULAR_MOMENTUM
298  )
299  );
300  return tabulated_real_quantities;
301  }
302 
303  void test_OrbitSolver::test_solution(
304  const std::vector< const std::list<double> * > &
305  tabulated_real_quantities,
306  std::vector<const Core::OneArgumentDiffFunction *>
307  expected_real_quantities,
308  const ExpectedEvolutionMode<Core::EvolModeType> &expected_evol_mode,
309  const ExpectedEvolutionMode<bool> &expected_wind_mode,
310  double min_age,
311  double max_age,
312  bool debug_mode
313  )
314  {
315  const std::list<Core::EvolModeType> &tabulated_modes =
316  __solver->mode_evolution();
317 
318  const std::list<bool> *tabulated_wind_sat =
319  (__star ? &(__star->wind_saturation_evolution()) : NULL);
320 
321  unsigned num_ages = tabulated_real_quantities[AGE]->size();
322 
323  std::ostringstream msg_start;
324  std::ostringstream msg;
325  msg.precision(16);
326  msg.setf(std::ios_base::scientific);
327  msg_start.precision(16);
328  msg_start.setf(std::ios_base::scientific);
329  msg << msg_start.str()
330  << num_ages
331  << " tabulated ages, ";
332  bool all_same_size = true;
333 
334  for(unsigned q = 0; q < AGE; ++q) {
335  if(tabulated_real_quantities[q]) {
336  msg << tabulated_real_quantities[q]->size()
337  << " tabulated "
338  << static_cast<RealEvolutionQuantity>(q)
339  << ", ";
340  all_same_size = (
341  all_same_size
342  &&
343  tabulated_real_quantities[q]->size() == num_ages
344  );
345  }
346  }
347 
348  msg << tabulated_modes.size() << " tabulated modes";
349 
350  all_same_size = (all_same_size
351  &&
352  tabulated_modes.size() == num_ages);
353  if(__star) {
354  msg << ", "
355  << tabulated_wind_sat->size()
356  << " tabulated wind saturations";
357  all_same_size = (all_same_size
358  &&
359  tabulated_wind_sat->size() == num_ages);
360  }
361 
362 
363  if(debug_mode) std::cout << msg.str() << std::endl;
364  TEST_ASSERT_MSG(all_same_size, msg.str().c_str());
365 
366  msg.str("");
367  msg << "Expected age range: " << min_age << " to " << max_age
368  << ", actual age range: " << tabulated_real_quantities[AGE]->front()
369  << " to " << tabulated_real_quantities[AGE]->back();
370  if(debug_mode) std::cout << msg.str() << std::endl;
371  TEST_ASSERT_MSG(
372  (
373  tabulated_real_quantities[AGE]->front() == min_age
374  &&
375  tabulated_real_quantities[AGE]->back() == max_age
376  ),
377  msg.str().c_str()
378  );
379 
380  std::vector< std::list<double>::const_iterator >
381  real_tabulated_iter(AGE);
382  for(unsigned q = 0; q < AGE; ++q)
383  if(tabulated_real_quantities[q])
384  real_tabulated_iter[q] = tabulated_real_quantities[q]->begin();
385  std::list<Core::EvolModeType>::const_iterator
386  tabulated_mode_iter = tabulated_modes.begin();
387  std::list<bool>::const_iterator tabulated_wind_sat_iter;
388  if(__star)
389  tabulated_wind_sat_iter = tabulated_wind_sat->begin();
390  double last_checked_age = Core::NaN;
391 
392  for(
393  std::list<double>::const_iterator
394  age_i = tabulated_real_quantities[AGE]->begin();
395  age_i != tabulated_real_quantities[AGE]->end();
396  ++age_i
397  ) {
398  std::vector<double> expected_real_values(AGE);
399  for(unsigned q = 0; q < AGE; ++q) {
400  expected_real_values[q] =
401  (*(expected_real_quantities[q]))(*age_i);
402  }
403  Core::EvolModeType expected_mode = expected_evol_mode(*age_i);
404  bool expected_wind_sat = expected_wind_mode(*age_i);
405 
406  std::ostringstream age_msg_start;
407  age_msg_start.precision(16);
408  age_msg_start.setf(std::ios_base::scientific);
409  age_msg_start << msg_start.str()
410  << "age = " << *age_i
411  << ", mode = " << *tabulated_mode_iter;
412  if(__star) {
413  age_msg_start << ", wind is ";
414  if(!(*tabulated_wind_sat_iter)) age_msg_start << " not ";
415  age_msg_start << "saturated";
416  }
417  for(unsigned q = 0; q < AGE; ++q)
418  if(tabulated_real_quantities[q])
419  age_msg_start << ", "
420  << static_cast<RealEvolutionQuantity>(q)
421  << " = "
422  << *real_tabulated_iter[q];
423 
424  msg.str("");
425  msg << age_msg_start.str() << " age is out of range.";
426  TEST_ASSERT_MSG(*age_i >= min_age && *age_i <= max_age,
427  msg.str().c_str());
428 
429  std::list<double>::const_iterator next_age_i = age_i;
430  ++next_age_i;
431  bool
432  can_skip = (
433  next_age_i != tabulated_real_quantities[AGE]->end()
434  &&
435  std::abs(*next_age_i - *age_i) < 1e-5
436  &&
437  expected_evol_mode.near_break(*age_i)
438  &&
439  expected_evol_mode.near_break(*next_age_i)
440  ),
441  skipped = (*age_i == last_checked_age);
442 
443  msg.str("");
444  msg << age_msg_start.str() << ": mode is not "
445  << expected_mode << ", but " << *tabulated_mode_iter;
446 
447  if(debug_mode) std::cout << msg.str() << std::endl;
448  if(!skipped && expected_mode != *tabulated_mode_iter) {
449  if(can_skip) skipped = true;
450  else TEST_ASSERT_MSG(false, msg.str().c_str());
451  }
452 
453  if(__star) {
454  msg.str("");
455  msg << age_msg_start.str() << ": wind is ";
456  if(!(*tabulated_wind_sat_iter)) msg << " not ";
457  msg << "saturated, but should";
458  if(!expected_wind_sat) msg << " not ";
459  msg << "be.";
460  if(debug_mode) std::cout << msg.str() << std::endl;
461  if(!skipped && expected_wind_sat != *tabulated_wind_sat_iter) {
462  if(can_skip) skipped = true;
463  else TEST_ASSERT_MSG(false, msg.str().c_str());
464  }
465  }
466 /* TEST_ASSERT_MSG(
467  (
468  skipped
469  ||
470  expected_wind_sat == *tabulated_wind_sat_iter
471  ||
472  expected_wind_mode.near_break(*age_i)
473  ),
474  msg.str().c_str()
475  );*/
476 
477  for(unsigned q = 0; q < AGE; ++q) {
478  if(!tabulated_real_quantities[q]) continue;
479  msg.str("");
480  msg << age_msg_start.str() << ": "
481  << static_cast<RealEvolutionQuantity>(q)
482  << " is not "
483  << expected_real_values[q]
484  << " but "
485  << *real_tabulated_iter[q]
486  << ", difference = "
487  << (*real_tabulated_iter[q]) - expected_real_values[q];
488  if(debug_mode) std::cout << msg.str() << std::endl;
489  if(
490  !skipped
491  &&
492  !check_diff((*real_tabulated_iter[q]),
493  expected_real_values[q],
494  1e-5,
495  0.0)
496  ) {
497  if(can_skip) {
498  skipped = true;
499  break;
500  }
501  TEST_ASSERT_MSG(false, msg.str().c_str());
502  }
503  }
504  if(!skipped)
505  last_checked_age = *age_i;
506  else if(debug_mode)
507  std::cerr << "Skipped checks for t = " << *age_i
508  << std::endl;
509  for(unsigned q = 0; q < AGE; ++q)
510  if(tabulated_real_quantities[q])
511  ++(real_tabulated_iter[q]);
512  if(__star) ++tabulated_wind_sat_iter;
513  ++tabulated_mode_iter;
514  }
515  }
516 
517  void test_OrbitSolver::test_no_planet_scenario(
518  const StellarEvolution::Interpolator &stellar_evol,
519  double *initial_Lstar,
520  double windK,
521  double wind_sat_freq,
522  double core_env_coupling_time,
523  std::vector<const Core::OneArgumentDiffFunction *>
524  &expected_real_quantities,
525  const ExpectedEvolutionMode<bool> &expected_wind_mode,
526  double max_age,
527  bool debug_mode
528  ) {
529  ExpectedEvolutionMode<Core::EvolModeType> single_mode,
530  binary_mode;
531  single_mode.add_break(TSTART, Core::SINGLE);
532  binary_mode.add_break(TSTART, Core::BINARY);
533 
534  __star = make_const_lag_star(stellar_evol,
535  windK,
536  wind_sat_freq,
537  core_env_coupling_time,
538  1.0);//phase lag
539  evolve(0.0,//Wdisk
540  0.0,//tdisk
541  1.0,//initial semimajor
542  initial_Lstar,
543  0.0,//initial inclination
544  1.0,//planet mass
545  max_age,//planet formation age.
546  max_age);//end evolution age
547  expected_real_quantities[SEMIMAJOR] = &nan_func;
548  expected_real_quantities[ECCENTRICITY] = &nan_func;
549  expected_real_quantities[CONV_INCLINATION] = &zero_func;
550  expected_real_quantities[RAD_INCLINATION] = &zero_func;
551  expected_real_quantities[CONV_PERIAPSIS] = &zero_func;
552  expected_real_quantities[RAD_PERIAPSIS] = &zero_func;
553 
554  test_solution(get_evolution(),
555  expected_real_quantities,
556  single_mode,
557  expected_wind_mode,
558  TSTART,
559  max_age,
560  debug_mode);
561 
562  if(initial_Lstar[0] == 0) return;
563 
564  for(double phase_lag = 0.0; phase_lag < 1.5; phase_lag += 1.0)
565  for(
566  double mplanet = 0.0;
567  mplanet < 1.5 - phase_lag;
568  mplanet += 1.0
569  ) {
570  delete __star;
571  delete __system;
572  delete __solver;
573 
574  __star = make_const_lag_star(stellar_evol,
575  windK,
576  wind_sat_freq,
577  core_env_coupling_time,
578  phase_lag);
579  evolve(0.0,//Wdisk
580  0.0,//tdisk
581  1.0,//initial semimajor
582  initial_Lstar,
583  0.0,//initial inclination
584  mplanet,//planet mass
585  TSTART,//planet formation age.
586  max_age);//end evolution age
587 
588  expected_real_quantities[SEMIMAJOR] = &one_func;
589  expected_real_quantities[ECCENTRICITY] = &zero_func;
590  test_solution(get_evolution(),
591  expected_real_quantities,
592  binary_mode,
593  expected_wind_mode,
594  TSTART,
595  max_age,
596  debug_mode);
597  delete __system;
598  delete __solver;
599 
600  evolve(0.0,//Wdisk
601  0.0,//tdisk
602  200.0,//initial semimajor
603  initial_Lstar,
604  0.0,//initial inclination
605  mplanet,//planet mass
606  TSTART,//planet formation age.
607  max_age);//end evolution age
608 
609  expected_real_quantities[SEMIMAJOR] = &two_hundred_func;
610  test_solution(get_evolution(),
611  expected_real_quantities,
612  binary_mode,
613  expected_wind_mode,
614  TSTART,
615  max_age,
616  debug_mode);
617 
618  }
619  delete __star;
620  __star = NULL;
621  delete __system;
622  __system = NULL;
623  delete __solver;
624  __solver = NULL;
625 
626  }
627 
628  std::vector<const Core::OneArgumentDiffFunction *>
629  test_OrbitSolver::calculate_expected_unlocked_evolution(
630  double phase_lag,
631  double secondary_mass,
632  bool decaying
633  )
634  {
635  std::vector<const Core::OneArgumentDiffFunction *>
636  expected_real_quantities(NUM_REAL_QUANTITIES - 1);
637 
638  double msecondary_si = (secondary_mass
639  *
640  Core::AstroConst::jupiter_mass),
641  mprimary_si = Core::AstroConst::solar_mass,
642  alpha = (
643  2.4 * M_PI
644  *
645  std::sqrt(
646  Core::AstroConst::G * (msecondary_si + mprimary_si)
647  /
648  Core::AstroConst::solar_radius
649  )
650  *
651  msecondary_si / mprimary_si
652  *
653  phase_lag * Core::AstroConst::Gyr
654  /
655  Core::AstroConst::solar_radius
656  ),
657  a6p5_offset,
658  a0,
659  Lscale = (
660  -msecondary_si
661  /
662  std::pow(Core::AstroConst::solar_radius, 1.5)
663  *
664  std::sqrt(
665  Core::AstroConst::G
666  /
667  (msecondary_si + Core::AstroConst::solar_mass)
668  )
669  *
670  Core::AstroConst::day
671  );
672 
673  if(decaying) {
674  a6p5_offset = std::pow(2.0, 6.5) + 6.5 * alpha;
675  a0 = std::pow(a6p5_offset, 1.0 / 6.5);
676  } else {
677  a0 = 2.6;
678  a6p5_offset = std::pow(a0, 6.5);
679  }
680 
681  std::valarray<double> a6p5_poly_coef(2);
682  a6p5_poly_coef[0] = a6p5_offset;
683  a6p5_poly_coef[1] = (decaying ? -1.0 : 1.0) * 6.5 * alpha;
684  StellarEvolution::PolynomialEvolutionQuantity *a6p5_evol = (
685  new StellarEvolution::PolynomialEvolutionQuantity(
686  a6p5_poly_coef,
687  TSTART,
688  1.0
689  )
690  );
691  __temp_functions.push_back(a6p5_evol);
692 
693  FunctionToPower *sqrta_evol = new FunctionToPower(a6p5_evol, 1.0/13.0);
694  __temp_functions.push_back(sqrta_evol);
695 
696 
697  ExponentialPlusFunc *Lconv_unscaled = new ExponentialPlusFunc(
698  sqrta_evol,
699  (decaying ? 0.0 : -1e5) - std::sqrt(a0),
700  0
701  );
702  __temp_functions.push_back(Lconv_unscaled);
703 
704  expected_real_quantities[SEMIMAJOR] = new FunctionToPower(
705  a6p5_evol,
706  1.0 / 6.5
707  );
708  __temp_functions.push_back(expected_real_quantities[SEMIMAJOR]);
709 
710  expected_real_quantities[ECCENTRICITY] = &zero_func;
711  expected_real_quantities[CONV_INCLINATION] = &zero_func;
712  expected_real_quantities[RAD_INCLINATION] = &zero_func;
713  expected_real_quantities[CONV_PERIAPSIS] = &zero_func;
714  expected_real_quantities[RAD_PERIAPSIS] = &zero_func;
715  expected_real_quantities[CONV_ANGMOM] = new ScaledFunction(
716  Lconv_unscaled,
717  Lscale
718  );
719  __temp_functions.push_back(expected_real_quantities[CONV_ANGMOM]);
720  expected_real_quantities[RAD_ANGMOM] = &zero_func;
721 
722  return expected_real_quantities;
723 
724  }
725 
726  std::vector<const Core::OneArgumentDiffFunction *>
727  test_OrbitSolver::calculate_expected_disklocked_to_fast_to_locked(
728  double lgQ,
729  double tdisk,
730  double async,
731  double tsync,
732  double tend,
733  bool include_disk_lock
734  )
735  {
736  const double Q = std::pow(10.0, lgQ),
737  alpha = (
738  -4.5
739  *
740  std::sqrt(
741  Core::AstroConst::G
742  /
743  (
744  Core::AstroConst::solar_radius
745  *
746  Core::AstroConst::solar_mass
747  )
748  )
749  *
750  Core::AstroConst::jupiter_mass / Q
751  *
752  Core::AstroConst::Gyr
753  /
754  Core::AstroConst::solar_radius
755  ),
756  Lscale = (
757  Core::AstroConst::jupiter_mass
758  /
759  std::pow(Core::AstroConst::solar_radius, 1.5)
760  *
761  std::sqrt(Core::AstroConst::G
762  /
763  (
764  Core::AstroConst::jupiter_mass
765  +
766  Core::AstroConst::solar_mass
767  )
768  )
769  *
770  Core::AstroConst::day
771  ),
772  beta = (
773  std::sqrt(
774  Core::AstroConst::G
775  *
776  (
777  Core::AstroConst::solar_mass
778  +
779  Core::AstroConst::jupiter_mass
780  )
781  )
782  *
783  Core::AstroConst::day
784  /
785  std::pow(Core::AstroConst::solar_radius, 1.5)
786  ),
787  a6p5_offset = (std::pow(async, 6.5)
788  -
789  6.5 * alpha * tsync),
790  a_formation=std::pow(
791  a6p5_offset + 6.5 * alpha * tdisk,
792  1.0 / 6.5
793  ),
794  Ic = (
795  Lscale
796  *
797  (std::sqrt(a_formation) - std::sqrt(async))
798  /
799  (beta * (std::pow(async, -1.5)
800  -
801  0.5 * std::pow(a_formation, -1.5)))
802  ),
803  wdisk = 0.5 * beta / std::pow(a_formation, 1.5),
804  wlocked = beta / std::pow(async, 1.5);
805 
806  std::valarray<double> a6p5_poly_coef(2);
807  a6p5_poly_coef[0] = a6p5_offset;
808  a6p5_poly_coef[1] = 6.5 * alpha;
809  StellarEvolution::PolynomialEvolutionQuantity
810  *a6p5_evol = new StellarEvolution::PolynomialEvolutionQuantity(
811  a6p5_poly_coef,
812  tdisk,
813  tsync
814  ),
815  *Lconv_disk = new StellarEvolution::PolynomialEvolutionQuantity(
816  std::valarray<double>(Ic * wdisk, 1),
817  TSTART,
818  tdisk
819  ),
820  *Lconv_locked = new StellarEvolution::PolynomialEvolutionQuantity(
821  std::valarray<double>(Ic * wlocked, 1),
822  tsync,
823  tend
824  ),
825  *nan_disk = new StellarEvolution::PolynomialEvolutionQuantity(
826  std::valarray<double>(Core::NaN, 2),
827  TSTART,
828  tdisk
829  ),
830  *a_locked = new StellarEvolution::PolynomialEvolutionQuantity(
831  std::valarray<double>(async, 1),
832  tsync,
833  tend
834  ),
835  *Lrad_evol = new StellarEvolution::PolynomialEvolutionQuantity(
836  std::valarray<double>(),
837  TSTART,
838  tend
839  ),
840  *zero_e = new StellarEvolution::PolynomialEvolutionQuantity(
841  std::valarray<double>(),
842  tdisk,
843  tend
844  );
845  __temp_functions.push_back(a6p5_evol);
846  __temp_functions.push_back(Lconv_disk);
847  __temp_functions.push_back(Lconv_locked);
848  __temp_functions.push_back(nan_disk);
849  __temp_functions.push_back(a_locked);
850  __temp_functions.push_back(Lrad_evol);
851  __temp_functions.push_back(zero_e);
852 
853  FunctionToPower *a_fast = new FunctionToPower(a6p5_evol,
854  1.0 / 6.5),
855  *sqrta_evol = new FunctionToPower(a6p5_evol,
856  1.0 / 13.0);
857  __temp_functions.push_back(a_fast);
858  __temp_functions.push_back(sqrta_evol);
859 
860  ExponentialPlusFunc *Lconv_unscaled = new ExponentialPlusFunc(
861  sqrta_evol,
862  -Ic * wdisk / Lscale - std::sqrt(a_formation), 0
863  );
864  __temp_functions.push_back(Lconv_unscaled);
865 
866  ScaledFunction *Lconv_fast = new ScaledFunction(Lconv_unscaled,
867  -Lscale);
868  __temp_functions.push_back(Lconv_fast);
869 
870  PiecewiseFunction *a_evol = new PiecewiseFunction,
871  *e_evol = new PiecewiseFunction,
872  *Lconv_evol = new PiecewiseFunction;
873  __temp_functions.push_back(a_evol);
874  __temp_functions.push_back(e_evol);
875  __temp_functions.push_back(Lconv_evol);
876 
877  StellarEvolution::PolynomialEvolutionQuantity *zero_quantity;
878 
879  if(include_disk_lock) {
880  a_evol->add_piece(nan_disk);
881  e_evol->add_piece(nan_disk);
882  Lconv_evol->add_piece(Lconv_disk);
883  zero_quantity = &zero_func;
884  } else {
885  zero_quantity = zero_e;
886  }
887 
888  a_evol->add_piece(a_fast);
889  a_evol->add_piece(a_locked);
890 
891  e_evol->add_piece(zero_e);
892 
893  Lconv_evol->add_piece(Lconv_fast);
894  Lconv_evol->add_piece(Lconv_locked);
895 
896  std::vector<const Core::OneArgumentDiffFunction *>
897  expected_real_quantities(NUM_REAL_QUANTITIES - 1);
898 
899  expected_real_quantities[SEMIMAJOR] = a_evol;
900  expected_real_quantities[ECCENTRICITY] = e_evol;
901  expected_real_quantities[CONV_INCLINATION] = zero_quantity;
902  expected_real_quantities[RAD_INCLINATION] = zero_quantity;
903  expected_real_quantities[CONV_PERIAPSIS] = zero_quantity;
904  expected_real_quantities[RAD_PERIAPSIS] = zero_quantity;
905  expected_real_quantities[CONV_ANGMOM] = Lconv_evol;
906  expected_real_quantities[RAD_ANGMOM] = zero_quantity;
907 
908  return expected_real_quantities;
909  }
910 
911  std::vector<const Core::OneArgumentDiffFunction *>
912  test_OrbitSolver::calculate_expected_polar_1_0(double tdisk,
913  double wstar,
914  double worb)
915 
916  {
917  double aorb = std::pow(
918  (
919  Core::AstroConst::G
920  *
921  (
922  Core::AstroConst::solar_mass
923  +
924  Core::AstroConst::jupiter_mass
925  )
926  )
927  /
928  std::pow(
929  worb / Core::AstroConst::day,
930  2
931  ),
932  1.0 / 3.0
933  ) / Core::AstroConst::solar_radius;
934 
935  StellarEvolution::PolynomialEvolutionQuantity
936  *disk_nan_evol = new StellarEvolution::PolynomialEvolutionQuantity(
937  std::valarray<double>(Core::NaN, 1),
938  TSTART,
939  tdisk
940  ),
941  *fixed_a_evol = new StellarEvolution::PolynomialEvolutionQuantity(
942  std::valarray<double>(aorb, 1),
943  tdisk,
944  MAX_AGE
945  ),
946  *fixed_e_evol = new StellarEvolution::PolynomialEvolutionQuantity(
947  std::valarray<double>(),
948  tdisk,
949  MAX_AGE
950  ),
951  *disk_zero_evol = new StellarEvolution::PolynomialEvolutionQuantity(
952  std::valarray<double>(),
953  TSTART,
954  tdisk
955  ),
956  *halfpi_evol = new StellarEvolution::PolynomialEvolutionQuantity(
957  std::valarray<double>(M_PI / 2.0, 1),
958  tdisk,
959  MAX_AGE
960  );
961  __temp_functions.push_back(disk_nan_evol);
962  __temp_functions.push_back(fixed_a_evol);
963  __temp_functions.push_back(fixed_e_evol);
964  __temp_functions.push_back(disk_zero_evol);
965  __temp_functions.push_back(halfpi_evol);
966 
967  PiecewiseFunction *a_evol = new PiecewiseFunction,
968  *e_evol = new PiecewiseFunction,
969  *conv_incl_evol = new PiecewiseFunction,
970  *rad_incl_evol = new PiecewiseFunction;
971  __temp_functions.push_back(a_evol);
972  __temp_functions.push_back(e_evol);
973  __temp_functions.push_back(conv_incl_evol);
974  __temp_functions.push_back(rad_incl_evol);
975 
976  a_evol->add_piece(disk_nan_evol);
977  a_evol->add_piece(fixed_a_evol);
978 
979  e_evol->add_piece(disk_nan_evol);
980  e_evol->add_piece(fixed_e_evol);
981 
982  conv_incl_evol->add_piece(disk_zero_evol);
983  conv_incl_evol->add_piece(halfpi_evol);
984 
985  rad_incl_evol->add_piece(disk_zero_evol);
986  rad_incl_evol->add_piece(halfpi_evol);
987 
988  std::vector<const Core::OneArgumentDiffFunction *>
989  expected_real_quantities(NUM_REAL_QUANTITIES - 1);
990 
991  expected_real_quantities[SEMIMAJOR] = a_evol;
992  expected_real_quantities[ECCENTRICITY] = e_evol;
993  expected_real_quantities[CONV_INCLINATION] = conv_incl_evol;
994  expected_real_quantities[RAD_INCLINATION] = rad_incl_evol;
995  expected_real_quantities[CONV_PERIAPSIS] = &zero_func;
996  expected_real_quantities[RAD_PERIAPSIS] = &zero_func;
997  expected_real_quantities[CONV_ANGMOM] =
998  new StellarEvolution::PolynomialEvolutionQuantity(
999  std::valarray<double>(wstar, 1),
1000  TSTART,
1001  MAX_AGE
1002  );
1003  __temp_functions.push_back(
1004  expected_real_quantities[CONV_ANGMOM]
1005  );
1006  expected_real_quantities[RAD_ANGMOM] = &one_func;
1007 
1008  return expected_real_quantities;
1009  }
1010 
1011  std::vector<const Core::OneArgumentDiffFunction *>
1012  test_OrbitSolver::calculate_expected_polar_2_0(double tdisk,
1013  double wstar,
1014  double worb,
1015  double phase_lag,
1016  double &lconv_decay_rate,
1017  double &semimajor)
1018  {
1019  semimajor = std::pow(
1020  (
1021  Core::AstroConst::G
1022  *
1023  (
1024  Core::AstroConst::solar_mass
1025  +
1026  Core::AstroConst::jupiter_mass
1027  )
1028  )
1029  /
1030  std::pow(
1031  worb / Core::AstroConst::day,
1032  2
1033  ),
1034  1.0 / 3.0
1035  ) / Core::AstroConst::solar_radius;
1036  lconv_decay_rate = (
1037  0.3 * M_PI
1038  *
1039  (
1040  Core::AstroConst::G
1041  *
1042  std::pow(
1043  (
1044  Core::AstroConst::jupiter_mass
1045  /
1046  std::pow(
1047  semimajor * Core::AstroConst::solar_radius,
1048  3
1049  )
1050  ),
1051  2
1052  )
1053  *
1054  std::pow(Core::AstroConst::solar_radius, 5)
1055  )
1056  /
1057  (
1058  Core::AstroConst::solar_mass
1059  *
1060  std::pow(Core::AstroConst::solar_radius, 2)
1061  )
1062  *
1063  Core::AstroConst::day
1064  *
1065  Core::AstroConst::Gyr
1066  *
1067  phase_lag
1068  );
1069  double lorb = (
1070  (
1071  Core::AstroConst::jupiter_mass
1072  /
1073  Core::AstroConst::solar_mass
1074  )
1075  /
1076  (
1077  1.0
1078  +
1079  Core::AstroConst::jupiter_mass
1080  /
1081  Core::AstroConst::solar_mass
1082  )
1083  *
1084  semimajor * semimajor
1085  *
1086  worb
1087  );
1088  double lconv_offset = lconv_decay_rate * tdisk;
1089 
1090  std::valarray<double> spindown_coef(2);
1091  spindown_coef[0] = wstar + lconv_decay_rate * tdisk;
1092  spindown_coef[1] = - lconv_decay_rate;
1093 
1094  std::valarray<double> conv_cosinc_numer_coef(3),
1095  conv_cosinc_denom_coef(2),
1096  rad_cosinc_numer_coef(3),
1097  rad_cosinc_denom_coef(3);
1098  conv_cosinc_numer_coef[0] = - (
1099  (2.0 * wstar + lconv_offset) * lconv_offset
1100  );
1101  conv_cosinc_numer_coef[1] = 2.0 * wstar * lconv_decay_rate;
1102  conv_cosinc_numer_coef[2] = -std::pow(lconv_decay_rate, 2);
1103 
1104  conv_cosinc_denom_coef[0] = 2.0 * lorb * (wstar + lconv_offset);
1105  conv_cosinc_denom_coef[1] = -2.0 * lconv_decay_rate * lorb;
1106 
1107  conv_cosinc_numer_coef[0] += conv_cosinc_denom_coef[0];
1108  conv_cosinc_numer_coef[1] += conv_cosinc_denom_coef[1];
1109 
1110  rad_cosinc_numer_coef[0] = - 2.0 * lorb * lconv_decay_rate * tdisk;
1111  rad_cosinc_numer_coef[1] = 2.0 * lorb * lconv_decay_rate;
1112 
1113  rad_cosinc_denom_coef[0] = (2.0 * lorb * lorb
1114  -
1115  2.0 * wstar * lconv_decay_rate * tdisk
1116  -
1117  std::pow(lconv_decay_rate * tdisk, 2));
1118  rad_cosinc_denom_coef[1] =
1119  2.0 * lconv_decay_rate * (worb + lconv_decay_rate * tdisk);
1120  rad_cosinc_denom_coef[2] = - std::pow(lconv_decay_rate, 2);
1121 
1122  rad_cosinc_numer_coef += rad_cosinc_denom_coef;
1123 
1124  StellarEvolution::PolynomialEvolutionQuantity
1125  *disk_nan_evol =
1126  new StellarEvolution::PolynomialEvolutionQuantity(
1127  std::valarray<double>(Core::NaN, 1),
1128  TSTART,
1129  tdisk
1130  ),
1131  *fixed_a_evol =
1132  new StellarEvolution::PolynomialEvolutionQuantity(
1133  std::valarray<double>(semimajor, 1),
1134  tdisk,
1135  MAX_AGE
1136  ),
1137  *fixed_e_evol
1138  = new StellarEvolution::PolynomialEvolutionQuantity(
1139  std::valarray<double>(),
1140  tdisk,
1141  MAX_AGE
1142  ),
1143  *disk_cosinc_evol =
1144  new StellarEvolution::PolynomialEvolutionQuantity(
1145  std::valarray<double>(2.0, 1),
1146  TSTART,
1147  tdisk
1148  ),
1149  *conv_cosinc_binary_numer =
1150  new StellarEvolution::PolynomialEvolutionQuantity(
1151  conv_cosinc_numer_coef,
1152  tdisk, MAX_AGE
1153  ),
1154  *conv_cosinc_binary_denom =
1155  new StellarEvolution::PolynomialEvolutionQuantity(
1156  conv_cosinc_denom_coef,
1157  tdisk,
1158  MAX_AGE
1159  ),
1160  *rad_cosinc_binary_numer =
1161  new StellarEvolution::PolynomialEvolutionQuantity(
1162  rad_cosinc_numer_coef,
1163  tdisk,
1164  MAX_AGE
1165  ),
1166  *rad_cosinc_binary_denom =
1167  new StellarEvolution::PolynomialEvolutionQuantity(
1168  rad_cosinc_denom_coef,
1169  tdisk,
1170  MAX_AGE
1171  ),
1172  *lconv_disk =
1173  new StellarEvolution::PolynomialEvolutionQuantity(
1174  std::valarray<double>(wstar, 1),
1175  TSTART,
1176  tdisk
1177  ),
1178  *lconv_spindown =
1179  new StellarEvolution::PolynomialEvolutionQuantity(
1180  spindown_coef,
1181  tdisk,
1182  MAX_AGE
1183  );
1184  __temp_functions.push_back(disk_nan_evol);
1185  __temp_functions.push_back(fixed_a_evol);
1186  __temp_functions.push_back(fixed_e_evol);
1187  __temp_functions.push_back(disk_cosinc_evol);
1188  __temp_functions.push_back(conv_cosinc_binary_numer);
1189  __temp_functions.push_back(conv_cosinc_binary_denom);
1190  __temp_functions.push_back(rad_cosinc_binary_numer);
1191  __temp_functions.push_back(rad_cosinc_binary_denom);
1192  __temp_functions.push_back(lconv_disk);
1193  __temp_functions.push_back(lconv_spindown);
1194 
1195  FunctionRatio *conv_cosinc_binary = new FunctionRatio(
1196  conv_cosinc_binary_numer,
1197  conv_cosinc_binary_denom
1198  );
1199  __temp_functions.push_back(conv_cosinc_binary);
1200  FunctionRatio *rad_cosinc_binary = new FunctionRatio(
1201  rad_cosinc_binary_numer,
1202  rad_cosinc_binary_denom
1203  );
1204  __temp_functions.push_back(rad_cosinc_binary);
1205 
1206  PiecewiseFunction *a_evol = new PiecewiseFunction,
1207  *e_evol = new PiecewiseFunction,
1208  *conv_cosinc_evol = new PiecewiseFunction,
1209  *rad_cosinc_evol = new PiecewiseFunction,
1210  *lconv_evol = new PiecewiseFunction;
1211  __temp_functions.push_back(a_evol);
1212  __temp_functions.push_back(e_evol);
1213  __temp_functions.push_back(conv_cosinc_evol);
1214  __temp_functions.push_back(rad_cosinc_evol);
1215  __temp_functions.push_back(lconv_evol);
1216 
1217  a_evol->add_piece(disk_nan_evol);
1218  a_evol->add_piece(fixed_a_evol);
1219 
1220  e_evol->add_piece(disk_nan_evol);
1221  e_evol->add_piece(fixed_e_evol);
1222 
1223  conv_cosinc_evol->add_piece(disk_cosinc_evol);
1224  conv_cosinc_evol->add_piece(conv_cosinc_binary);
1225 
1226  rad_cosinc_evol->add_piece(disk_cosinc_evol);
1227  rad_cosinc_evol->add_piece(rad_cosinc_binary);
1228 
1229  lconv_evol->add_piece(lconv_disk);
1230  lconv_evol->add_piece(lconv_spindown);
1231 
1232  std::vector<const Core::OneArgumentDiffFunction *>
1233  expected_real_quantities(NUM_REAL_QUANTITIES - 1);
1234 
1235  expected_real_quantities[SEMIMAJOR] = a_evol;
1236  expected_real_quantities[ECCENTRICITY] = e_evol;
1237  expected_real_quantities[CONV_INCLINATION] = conv_cosinc_evol;
1238  expected_real_quantities[RAD_INCLINATION] = rad_cosinc_evol;
1239  expected_real_quantities[CONV_PERIAPSIS] = &zero_func;
1240  expected_real_quantities[RAD_PERIAPSIS] = &zero_func;
1241  expected_real_quantities[CONV_ANGMOM] = lconv_evol;
1242  expected_real_quantities[RAD_ANGMOM] = &one_func;
1243 
1244  return expected_real_quantities;
1245  }
1246 
1247  std::vector<const Core::OneArgumentDiffFunction *>
1248  test_OrbitSolver::calculate_expected_oblique_m_0(
1249  unsigned m,
1250  double tdisk,
1251  double worb,
1252  double initial_inc,
1253  double initial_wstar,
1254  double phase_lag,
1255  double &min_wstar
1256  )
1257  {
1258  assert(m == 1 || m == 2);
1259  double aorb = std::pow(
1260  (
1261  Core::AstroConst::G
1262  *
1263  (
1264  Core::AstroConst::solar_mass
1265  +
1266  Core::AstroConst::jupiter_mass
1267  )
1268  )
1269  /
1270  std::pow(
1271  worb / Core::AstroConst::day,
1272  2
1273  ),
1274  1.0 / 3.0
1275  ) / Core::AstroConst::solar_radius;
1276  double lorb = (Mjup_to_Msun
1277  /
1278  (1.0 + Mjup_to_Msun)
1279  *
1280  aorb * aorb
1281  *
1282  worb);
1283 
1284  double ltot = std::sqrt(
1285  std::pow(initial_wstar, 2)
1286  +
1287  2.0 * initial_wstar * lorb * std::cos(initial_inc)
1288  +
1289  std::pow(lorb, 2)
1290  );
1291 
1292  min_wstar = ltot - lorb;
1293 
1294  double linear_quantity_rate = (
1295  0.6 * M_PI
1296  *
1297  (
1298  Core::AstroConst::G
1299  *
1300  std::pow(
1301  (
1302  Core::AstroConst::jupiter_mass
1303  /
1304  std::pow(
1305  aorb * Core::AstroConst::solar_radius,
1306  3
1307  )
1308  ),
1309  2
1310  )
1311  *
1312  std::pow(Core::AstroConst::solar_radius, 5)
1313  )
1314  /
1315  (
1316  Core::AstroConst::solar_mass
1317  *
1318  std::pow(Core::AstroConst::solar_radius, 2)
1319  )
1320  *
1321  Core::AstroConst::day
1322  *
1323  Core::AstroConst::Gyr
1324  *
1325  phase_lag
1326  );
1327 
1328  StellarEvolution::PolynomialEvolutionQuantity *disk_nan_evol =
1329  new StellarEvolution::PolynomialEvolutionQuantity(
1330  std::valarray<double>(Core::NaN, 1),
1331  TSTART,
1332  tdisk
1333  );
1334  __temp_functions.push_back(disk_nan_evol);
1335 
1336  StellarEvolution::PolynomialEvolutionQuantity *fixed_a_evol =
1337  new StellarEvolution::PolynomialEvolutionQuantity(
1338  std::valarray<double>(aorb, 1),
1339  tdisk,
1340  MAX_AGE
1341  );
1342  __temp_functions.push_back(fixed_a_evol);
1343 
1344  StellarEvolution::PolynomialEvolutionQuantity *fixed_e_evol =
1345  new StellarEvolution::PolynomialEvolutionQuantity(
1346  std::valarray<double>(),
1347  tdisk,
1348  MAX_AGE
1349  );
1350  __temp_functions.push_back(fixed_e_evol);
1351 
1352  PiecewiseFunction *a_evol = new PiecewiseFunction,
1353  *e_evol = new PiecewiseFunction;
1354  __temp_functions.push_back(a_evol);
1355  __temp_functions.push_back(e_evol);
1356 
1357  a_evol->add_piece(disk_nan_evol);
1358  a_evol->add_piece(fixed_a_evol);
1359 
1360  e_evol->add_piece(disk_nan_evol);
1361  e_evol->add_piece(fixed_e_evol);
1362 
1363  Core::OneArgumentDiffFunction *lconv_evol;
1364  if(m == 1)
1365  lconv_evol =
1366  new InverseLinearLconvEvolution<Oblique10LinearQuantity>(
1367  tdisk,
1368  ltot,
1369  lorb,
1370  initial_wstar,
1371  linear_quantity_rate
1372  );
1373  else
1374  lconv_evol =
1375  new InverseLinearLconvEvolution<Oblique20LinearQuantity>(
1376  tdisk,
1377  ltot,
1378  lorb,
1379  initial_wstar,
1380  linear_quantity_rate / 2.0
1381  );
1382  __temp_functions.push_back(lconv_evol);
1383 
1384  ConservedLEConvObliquityEvolution *conv_obliq_evol =
1385  new ConservedLEConvObliquityEvolution(*lconv_evol,
1386  lorb,
1387  ltot,
1388  tdisk);
1389  __temp_functions.push_back(conv_obliq_evol);
1390 
1391  ConservedLERadObliquityEvolution *rad_obliq_evol =
1392  new ConservedLERadObliquityEvolution(*lconv_evol,
1393  *conv_obliq_evol,
1394  lorb,
1395  tdisk);
1396  __temp_functions.push_back(rad_obliq_evol);
1397 
1398  std::vector<const Core::OneArgumentDiffFunction *>
1399  expected_real_quantities(NUM_REAL_QUANTITIES - 1);
1400  expected_real_quantities[SEMIMAJOR] = a_evol;
1401  expected_real_quantities[ECCENTRICITY] = e_evol;
1402  expected_real_quantities[CONV_INCLINATION] = conv_obliq_evol;
1403  expected_real_quantities[RAD_INCLINATION] = rad_obliq_evol;
1404  expected_real_quantities[CONV_PERIAPSIS] = &zero_func;
1405  expected_real_quantities[RAD_PERIAPSIS] = &zero_func;
1406  expected_real_quantities[CONV_ANGMOM] = lconv_evol;
1407  expected_real_quantities[RAD_ANGMOM] = &one_func;
1408 
1409  return expected_real_quantities;
1410  }
1411 
1412  void test_OrbitSolver::test_disk_locked_no_stellar_evolution()
1413  {
1414  try {
1415  StellarEvolution::MockStellarEvolution *no_evol =
1416  StellarEvolution::make_no_evolution();
1417  __star = make_const_lag_star(*no_evol,
1418  1.0,
1419  1.0,
1420  1.0);
1421 
1422  ExpectedEvolutionMode<Core::EvolModeType> expected_evol_mode;
1423  expected_evol_mode.add_break(TSTART, Core::LOCKED_SURFACE_SPIN);
1424 
1425  ExpectedEvolutionMode<bool> expected_wind_mode;
1426  expected_wind_mode.add_break(TSTART, false);
1427 
1428  std::vector<const Core::OneArgumentDiffFunction *>
1429  expected_real_quantities(NUM_REAL_QUANTITIES - 1);
1430  expected_real_quantities[SEMIMAJOR] = &nan_func;
1431  expected_real_quantities[ECCENTRICITY] = &nan_func;
1432  expected_real_quantities[CONV_INCLINATION] = &zero_func;
1433  expected_real_quantities[RAD_INCLINATION] = &zero_func;
1434  expected_real_quantities[CONV_PERIAPSIS] = &zero_func;
1435  expected_real_quantities[RAD_PERIAPSIS] = &zero_func;
1436  expected_real_quantities[CONV_ANGMOM] = &one_func;
1437  expected_real_quantities[RAD_ANGMOM] = &one_func;
1438 
1439  evolve(1.0, MAX_AGE, 1.0, &one);
1440  test_solution(get_evolution(),
1441  expected_real_quantities,
1442  expected_evol_mode,
1443  expected_wind_mode,
1444  TSTART,
1445  MAX_AGE);
1446 
1447  delete __solver;
1448  delete __system;
1449 
1450  evolve(1.0, MAX_AGE, 1.0, &zero);
1451  expected_real_quantities[RAD_ANGMOM] = new ExponentialPlusFunc(
1452  &one_func,
1453  -std::exp(TSTART/2.0),
1454  -0.5
1455  );
1456  test_solution(get_evolution(),
1457  expected_real_quantities,
1458  expected_evol_mode,
1459  expected_wind_mode,
1460  TSTART,
1461  MAX_AGE);
1462 
1463  delete expected_real_quantities[RAD_ANGMOM];
1464  delete no_evol;
1465  } catch (Core::Error::General &ex) {
1466  TEST_ASSERT_MSG(false, (std::string("Unexpected exception thrown: ")+
1467  ex.what()+": "+ex.get_message()).c_str());
1468  } catch (std::exception &ex) {
1469  TEST_ASSERT_MSG(false, (std::string("Unexpected exception thrown: ")+
1470  ex.what()).c_str());
1471  }
1472  }
1473 
1474  void test_OrbitSolver::test_disk_locked_with_stellar_evolution()
1475  {
1476  try {
1477  StellarEvolution::MockStellarEvolution *evol1 =
1478  StellarEvolution::make_linear_I_evolution();
1479 
1480  __star = make_const_lag_star(*evol1, 1.0, 1.0, 1.0);
1481 
1482  ExpectedEvolutionMode<Core::EvolModeType> expected_evol_mode;
1483  expected_evol_mode.add_break(TSTART, Core::LOCKED_SURFACE_SPIN);
1484 
1485  ExpectedEvolutionMode<bool> expected_wind_mode;
1486  expected_wind_mode.add_break(TSTART, false);
1487 
1488  std::valarray<double> temp_array=std::valarray<double>(1.0, 2);
1489  temp_array[0] = -1;
1490  StellarEvolution::PolynomialEvolutionQuantity *temp_poly =
1491  new StellarEvolution::PolynomialEvolutionQuantity(temp_array,
1492  TSTART,
1493  MAX_AGE);
1494 
1495  std::vector<const Core::OneArgumentDiffFunction *>
1496  expected_real_quantities(NUM_REAL_QUANTITIES - 1);
1497  expected_real_quantities[SEMIMAJOR] = &nan_func;
1498  expected_real_quantities[ECCENTRICITY] = &nan_func;
1499  expected_real_quantities[CONV_INCLINATION] = &zero_func;
1500  expected_real_quantities[RAD_INCLINATION] = &zero_func;
1501  expected_real_quantities[CONV_PERIAPSIS] = &zero_func;
1502  expected_real_quantities[RAD_PERIAPSIS] = &zero_func;
1503  expected_real_quantities[CONV_ANGMOM] =
1504  new StellarEvolution::PolynomialEvolutionQuantity(
1505  std::valarray<double>(1.0, 2),
1506  TSTART,
1507  MAX_AGE
1508  );
1509  expected_real_quantities[RAD_ANGMOM] =
1510  new ExponentialPlusFunc(temp_poly, 1, -0.5);
1511 
1512  double initial_Lrad = TSTART - 1 + std::exp(-TSTART / 2.0);
1513  evolve(1.0,
1514  MAX_AGE,
1515  1.0,
1516  &initial_Lrad);
1517  test_solution(get_evolution(),
1518  expected_real_quantities,
1519  expected_evol_mode,
1520  expected_wind_mode,
1521  TSTART,
1522  MAX_AGE);
1523 
1524  delete expected_real_quantities[CONV_ANGMOM];
1525  delete expected_real_quantities[RAD_ANGMOM];
1526  delete temp_poly;
1527  delete __solver;
1528  delete __system;
1529  delete __star;
1530  delete evol1;
1531 
1532  std::valarray< std::valarray<double> > Ic_arr(
1533  std::valarray<double>(1.0, 1),
1534  2
1535  );
1536  Ic_arr[1][0]=-1.0/6.0;
1537  StellarEvolution::MockStellarEvolution evol2(
1538  -1,
1539  std::valarray< std::valarray<double> >(//R
1540  std::valarray<double>(1.0, 1),
1541  1
1542  ),
1543  Ic_arr, //Iconv
1544  std::valarray< std::valarray<double> >(//Irad
1545  std::valarray<double>(1.0/6.0, 1),
1546  2
1547  ),
1548  std::valarray< std::valarray<double> >(//Rcore
1549  std::valarray<double>(0.5, 1),
1550  1
1551  ),
1552  std::valarray< std::valarray<double> >(//Mcore
1553  std::valarray<double>(1.0, 1),
1554  2
1555  )
1556  );
1557 
1558  __star = make_const_lag_star(evol2, 1.0, 1.0, 1.0);
1559 
1560  std::valarray<double> Lc_coef(1.0, 2);
1561  Lc_coef[1]=-1.0/6.0;
1562  expected_real_quantities[SEMIMAJOR] = &nan_func;
1563  expected_real_quantities[ECCENTRICITY] = &nan_func;
1564  expected_real_quantities[CONV_INCLINATION] = &zero_func;
1565  expected_real_quantities[RAD_INCLINATION] = &zero_func;
1566  expected_real_quantities[CONV_PERIAPSIS] = &zero_func;
1567  expected_real_quantities[RAD_PERIAPSIS] = &zero_func;
1568  expected_real_quantities[CONV_ANGMOM] =
1569  new StellarEvolution::PolynomialEvolutionQuantity(Lc_coef,
1570  TSTART,
1571  MAX_AGE);
1572  expected_real_quantities[RAD_ANGMOM] =
1573  new StellarEvolution::PolynomialEvolutionQuantity(
1574  std::valarray<double>(1.0/6.0, 2),
1575  TSTART,
1576  MAX_AGE
1577  );
1578  initial_Lrad = (TSTART / 2.0 + 1.0) / 6.0;
1579  evolve(1.0,
1580  MAX_AGE,
1581  1.0,
1582  &initial_Lrad);
1583  test_solution(get_evolution(),
1584  expected_real_quantities,
1585  expected_evol_mode,
1586  expected_wind_mode,
1587  TSTART,
1588  MAX_AGE);
1589 
1590  delete expected_real_quantities[CONV_ANGMOM];
1591  delete expected_real_quantities[RAD_ANGMOM];
1592  } catch (Core::Error::General &ex) {
1593  TEST_ASSERT_MSG(false, (std::string("Unexpected exception thrown: ")
1594  +
1595  ex.what()+": "+ex.get_message()).c_str());
1596  } catch (std::exception &ex) {
1597  TEST_ASSERT_MSG(false, (std::string("Unexpected exception thrown: ")
1598  +
1599  ex.what()).c_str());
1600  }
1601  }
1602 
1603  void test_OrbitSolver::test_no_planet_evolution()
1604  {
1605  try {
1606  const double rt2 = std::sqrt(2.0);
1607 
1608  StellarEvolution::MockStellarEvolution
1609  *stellar_evol = StellarEvolution::make_no_evolution();
1610  double initial_Lstar[] = {0.0, 0.0};
1611 
1612  std::vector<const Core::OneArgumentDiffFunction *>
1613  expected_real_quantities(NUM_REAL_QUANTITIES - 1);
1614  expected_real_quantities[CONV_ANGMOM] = &zero_func;
1615  expected_real_quantities[RAD_ANGMOM] = &zero_func;
1616  ExpectedEvolutionMode<bool> unsat_wind_mode, sat_wind_mode;
1617  unsat_wind_mode.add_break(TSTART, false);
1618  sat_wind_mode.add_break(TSTART, true);
1619 
1620  test_no_planet_scenario(*stellar_evol,
1621  initial_Lstar,
1622  1.0,//Wind K
1623  2.0,//wind sat freq.
1624  Core::Inf,//core-env coupling timescale
1625  expected_real_quantities,
1626  unsat_wind_mode);
1627 
1628  delete stellar_evol;
1629  stellar_evol = StellarEvolution::make_linear_I_evolution();
1630 
1631  initial_Lstar[0] = 1.0;
1632  expected_real_quantities[CONV_ANGMOM] = &one_func;
1633  test_no_planet_scenario(*stellar_evol,
1634  initial_Lstar,
1635  0.0,//Wind K
1636  2.0,//wind sat freq.
1637  Core::Inf,//core-env coupling timescale
1638  expected_real_quantities,
1639  unsat_wind_mode);
1640 
1641 
1642  initial_Lstar[0] = 0.5 * (1.0 + std::exp(-TSTART));
1643  initial_Lstar[1] = 0.5 * (1.0 - std::exp(-TSTART));
1644 
1645  StellarEvolution::PolynomialEvolutionQuantity half_func(
1646  std::valarray<double>(0.5, 1),
1647  TSTART,
1648  MAX_AGE
1649  );
1650  expected_real_quantities[CONV_ANGMOM] = new ExponentialPlusFunc(
1651  &half_func,
1652  0.5,
1653  -1.0
1654  );
1655  expected_real_quantities[RAD_ANGMOM] = new ExponentialPlusFunc(
1656  &half_func,
1657  -0.5,
1658  -1.0
1659  );
1660  test_no_planet_scenario(*stellar_evol,
1661  initial_Lstar,
1662  0.0,//Wind K
1663  1.0,//wind sat freq.
1664  1.0,//core-env coupling timescale
1665  expected_real_quantities,
1666  unsat_wind_mode);
1667 
1668  delete expected_real_quantities[CONV_ANGMOM];
1669  delete expected_real_quantities[RAD_ANGMOM];
1670 
1671  initial_Lstar[0] = 1.0/(1.0+TSTART);
1672  initial_Lstar[1] = 1.0;
1673 
1674  double wind_sat_age = std::pow(2.0, 0.25) - 1;
1675  std::valarray<double> late_denom_coef(3);
1676  late_denom_coef[0] = 2.0 * rt2 - 2.0;
1677  late_denom_coef[1] = 4.0 * rt2;
1678  late_denom_coef[2] = 2.0 * rt2;
1679  StellarEvolution::PolynomialEvolutionQuantity
1680  one_func_early(std::valarray<double>(1.0, 1),
1681  TSTART, wind_sat_age),
1682  one_plus_t(std::valarray<double>(1.0, 2), TSTART, MAX_AGE),
1683  late_denom2(late_denom_coef, wind_sat_age, MAX_AGE);
1684  FunctionToPower late_denom(&late_denom2, 0.5);
1685  FunctionRatio early_solution(&one_func_early, &one_plus_t),
1686  late_solution(&one_plus_t, &late_denom);
1687  PiecewiseFunction full_solution;
1688  full_solution.add_piece(&early_solution);
1689  full_solution.add_piece(&late_solution);
1690 
1691  expected_real_quantities[CONV_ANGMOM] = &full_solution;
1692  expected_real_quantities[RAD_ANGMOM] = &one_func;
1693  ExpectedEvolutionMode<bool> changing_wind_mode;
1694  changing_wind_mode.add_break(TSTART, true);
1695  changing_wind_mode.add_break(wind_sat_age, false);
1696 
1697  test_no_planet_scenario(*stellar_evol,
1698  initial_Lstar,
1699  2.0,//Wind K
1700  1.0 / rt2, //wind sat freq.
1701  Core::Inf,//core-env coupling timescale
1702  expected_real_quantities,
1703  changing_wind_mode);
1704 
1705  double b1 = (std::sqrt(2) - 1) / (2.0 * rt2),
1706  b2 = (std::sqrt(2) + 1) / (2.0 * rt2);
1707  initial_Lstar[0] = (b1 * std::exp(-TSTART / (2.0 + rt2))
1708  +
1709  b2 * std::exp(-TSTART / (2.0 - rt2)));
1710  initial_Lstar[1] = (
1711  0.5 / rt2 * std::exp(-1.0 / (2.0 + rt2) * (TSTART))
1712  -
1713  0.5 / rt2 * std::exp(-1.0 / (2.0 - rt2) * (TSTART))
1714  );
1715 
1716  ExponentialPlusFunc
1717  Lc1(&zero_func, b1, -1.0 / (2.0 + rt2)),
1718  Lc2(&zero_func, b2, -1.0 / (2.0 - rt2)),
1719  Lr1(&zero_func, 0.5 / rt2, -1.0 / (2.0 + rt2)),
1720  Lr2(&zero_func, -0.5 / rt2, -1.0 / (2.0 - rt2));
1721 
1722  expected_real_quantities[CONV_ANGMOM] = new FuncPlusFunc(&Lc1,
1723  &Lc2);
1724  expected_real_quantities[RAD_ANGMOM] = new FuncPlusFunc(&Lr1,
1725  &Lr2);
1726  delete stellar_evol;
1727  stellar_evol = StellarEvolution::make_no_evolution();
1728  test_no_planet_scenario(*stellar_evol,
1729  initial_Lstar,
1730  100.0,//Wind K
1731  0.1,//wind sat freq.
1732  1.0,//core-env coupling timescale
1733  expected_real_quantities,
1734  sat_wind_mode,
1735  2.0);
1736  delete expected_real_quantities[CONV_ANGMOM];
1737  delete expected_real_quantities[RAD_ANGMOM];
1738 
1739  delete stellar_evol;
1740  } catch (Core::Error::General &ex) {
1741  TEST_ASSERT_MSG(false, (std::string("Unexpected exception thrown: ")+
1742  ex.what()+": "+ex.get_message()).c_str());
1743  } catch (std::exception &ex) {
1744  TEST_ASSERT_MSG(false, (std::string("Unexpected exception thrown: ")+
1745  ex.what()).c_str());
1746  }
1747  }
1748 
1749  void test_OrbitSolver::test_unlocked_evolution()
1750  {
1751  try {
1752  ExpectedEvolutionMode<Core::EvolModeType> binary_mode;
1753  binary_mode.add_break(TSTART, Core::BINARY);
1754  ExpectedEvolutionMode<bool> unsat_wind_mode, sat_wind_mode;
1755  unsat_wind_mode.add_break(TSTART, false);
1756  sat_wind_mode.add_break(TSTART, true);
1757 
1758  StellarEvolution::MockStellarEvolution
1759  *no_evol = StellarEvolution::make_no_evolution(1.0);
1760 
1761  const double mplanet = 100;
1762  double lag = 1e-8 / mplanet;
1763 
1764  std::vector<const Core::OneArgumentDiffFunction *>
1765  expected_real_quantities = calculate_expected_unlocked_evolution(
1766  lag,
1767  mplanet
1768  );
1769 
1770  __star = make_const_lag_star(*no_evol,
1771  0.0,//wind K
1772  100.0,//wsat
1773  Core::Inf,//tcoup
1774  lag);
1775 
1776  while(true) {
1777  double initial_a = (
1778  *expected_real_quantities[SEMIMAJOR]
1779  )(TSTART);
1780 
1781  std::valarray<double> initial_L(3);
1782  initial_L[0] = (*expected_real_quantities[CONV_ANGMOM])(
1783  TSTART
1784  );
1785  initial_L[1] = 0.0;
1786  initial_L[2] = 0.0;
1787 
1788  evolve(0.0,//wdisk
1789  TSTART,//tdisk
1790  initial_a,
1791  &(initial_L[0]),
1792  0.0,//initial inclination
1793  mplanet,//planet mass
1794  Core::NaN,//tplanet
1795  1.0,//stop evolution age
1796  0.0001);//Rplanet
1797  test_solution(get_evolution(),
1798  expected_real_quantities,
1799  binary_mode,
1800  unsat_wind_mode,
1801  TSTART,
1802  1.0);
1803 
1804  delete __system;
1805  delete __solver;
1806 
1807  expected_real_quantities =
1808  calculate_expected_unlocked_evolution(
1809  lag,
1810  mplanet,
1811  false
1812  );
1813 
1814  initial_a = (*expected_real_quantities[SEMIMAJOR])(TSTART);
1815  initial_L[0] = (*expected_real_quantities[CONV_ANGMOM])(TSTART);
1816 
1817  evolve(0.0,//wdisk
1818  TSTART,//tdisk
1819  initial_a,
1820  &(initial_L[0]),
1821  0.0,//initial inclination
1822  mplanet,//planet mass
1823  Core::NaN,//tplanet
1824  1.0,//stop evolution age
1825  0.0001);//Rplanet
1826  test_solution(get_evolution(),
1827  expected_real_quantities,
1828  binary_mode,
1829  sat_wind_mode,
1830  TSTART,
1831  1.0);
1832 
1833  if(__star) {
1834  delete __star;
1835  __star = NULL;
1836 
1837  __primary_planet = new Planet::Planet(1.0, 1.0, 1.0);
1838  __primary_planet->zone().setup(
1839  std::vector<double>(),//Wtide breaks
1840  std::vector<double>(),//W* breaks
1841  std::vector<double>(1, 0.0),//Wtide pow.
1842  std::vector<double>(1, 0.0),//W* pow.
1843  lag
1844  );
1845  } else
1846  break;
1847  }
1848 
1849  delete no_evol;
1850  } catch (Core::Error::General &ex) {
1851  TEST_ASSERT_MSG(false, (std::string("Unexpected exception thrown: ")
1852  +
1853  ex.what()+": "+ex.get_message()).c_str());
1854  } catch (std::exception &ex) {
1855  TEST_ASSERT_MSG(false, (std::string("Unexpected exception thrown: ")
1856  +
1857  ex.what()).c_str());
1858  }
1859  }
1860 
1865  double locked_unsat_eq(double a, void *params)
1866  {
1867  double *dbl_par = static_cast<double *>(params);
1868  double alpha = dbl_par[0],
1869  beta = dbl_par[1],
1870  kappa = dbl_par[2],
1871  c = dbl_par[3],
1872  t = dbl_par[4];
1873  return (
1874  0.1 * alpha * std::pow(a, 5)
1875  -
1876  0.5 * beta * std::pow(a, 3)
1877  +
1878  kappa*t
1879  -
1880  c
1881  );
1882  }
1883 
1888  double locked_sat_eq(double a, void *params)
1889  {
1890  double *dbl_par = static_cast<double *>(params);
1891  double alpha = dbl_par[0],
1892  beta = dbl_par[1],
1893  kappa = dbl_par[2],
1894  c = dbl_par[3],
1895  t = dbl_par[4];
1896  return (
1897  alpha * a * a / 4.0
1898  -
1899  1.5 * beta * std::log(a)
1900  +
1901  kappa * t
1902  -
1903  c
1904  );
1905  }
1906 
1912  double locked_unsat_deriv(double a, void *params)
1913  {
1914  double *dbl_par = static_cast<double *>(params);
1915  double alpha = dbl_par[0],
1916  beta = dbl_par[1];
1917  return 0.5 * alpha * std::pow(a, 4) - 1.5 * beta * std::pow(a, 2);
1918  }
1919 
1925  double locked_sat_deriv(double a, void *params)
1926  {
1927  double *dbl_par = static_cast<double *>(params);
1928  double alpha = dbl_par[0],
1929  beta = dbl_par[1];
1930  return alpha * a / 2.0 - 1.5 * beta / a;
1931  }
1932 
1938  void locked_unsat_eq_deriv(double a, void *params, double *f, double *df)
1939  {
1940  double *dbl_par = static_cast<double *>(params);
1941  double alpha = dbl_par[0],
1942  beta = dbl_par[1],
1943  kappa = dbl_par[2],
1944  c = dbl_par[3],
1945  t = dbl_par[4];
1946  *f = (
1947  0.1 * alpha * std::pow(a, 5)
1948  -
1949  0.5 * beta * std::pow(a, 3)
1950  +
1951  kappa * t
1952  -
1953  c
1954  );
1955  *df = (0.5 * alpha * std::pow(a, 4)
1956  -
1957  1.5 * beta * std::pow(a, 2));
1958  }
1959 
1965  void locked_sat_eq_deriv(double a, void *params, double *f, double *df)
1966  {
1967  double *dbl_par = static_cast<double *>(params);
1968  double alpha = dbl_par[0],
1969  beta = dbl_par[1],
1970  kappa = dbl_par[2],
1971  c = dbl_par[3],
1972  t = dbl_par[4];
1973  *f = (
1974  alpha * a * a / 4.0
1975  -
1976  1.5 * beta * std::log(a)
1977  +
1978  kappa * t
1979  -
1980  c
1981  );
1982  *df = alpha * a / 2.0 - 1.5 * beta / a;
1983  }
1984 
1985  void test_OrbitSolver::test_locked_evolution()
1986  {
1987  try {
1988  ExpectedEvolutionMode<Core::EvolModeType> expected_mode;
1989  expected_mode.add_break(TSTART, Core::BINARY);
1990  ExpectedEvolutionMode<bool> expected_wind_mode;
1991  expected_wind_mode.add_break(TSTART, false);
1992  const double Ic = 0.001, Kwind = 1e-3, Kwind_s = 1;
1993  StellarEvolution::MockStellarEvolution *
1994  no_evol = StellarEvolution::make_no_evolution(1.0, Ic);
1995 
1996  double a1 = 3,
1997  Lscale = (Core::AstroConst::jupiter_mass
1998  /
1999  std::pow(Core::AstroConst::solar_radius, 1.5)
2000  *
2001  std::sqrt(
2002  Core::AstroConst::G
2003  /
2004  (Core::AstroConst::jupiter_mass
2005  +
2006  Core::AstroConst::solar_mass)
2007  )
2008  *
2009  Core::AstroConst::day),
2010  beta = (
2011  Ic * std::sqrt(
2012  Core::AstroConst::G
2013  *
2014  (
2015  Core::AstroConst::solar_mass
2016  +
2017  Core::AstroConst::jupiter_mass
2018  )
2019  )
2020  *
2021  Core::AstroConst::day
2022  /
2023  std::pow(Core::AstroConst::solar_radius, 1.5)
2024  ),
2025  wsat_s = 0.1, //Must be adjusted if a1 is adjusted
2026 // wsat = Core::Inf,
2027  kappa=(
2028  Kwind
2029  *
2030  std::pow(
2031  (
2032  Core::AstroConst::G
2033  *
2034  (
2035  Core::AstroConst::solar_mass
2036  +
2037  Core::AstroConst::jupiter_mass
2038  )
2039  /
2040  std::pow(Core::AstroConst::solar_radius, 3)
2041  ),
2042  1.5
2043  )
2044  *
2045  std::pow(Core::AstroConst::day, 3)
2046  ),
2047  kappa_s = (
2048  Kwind_s * wsat_s * wsat_s *
2049  std::sqrt(
2050  Core::AstroConst::G
2051  *
2052  (
2053  Core::AstroConst::solar_mass
2054  +
2055  Core::AstroConst::jupiter_mass
2056  )
2057  /
2058  std::pow(Core::AstroConst::solar_radius, 3)
2059  )
2060  *
2061  Core::AstroConst::day
2062  ),
2063  int_const = (Lscale / 10.0 * std::pow(a1, 5)
2064  -
2065  beta / 2.0 * std::pow(a1, 3)
2066  +
2067  kappa),
2068  int_const_s = (Lscale / 4.0 * std::pow(a1, 2)
2069  -
2070  3.0 * beta / 2.0 * std::log(a1)
2071  +
2072  kappa_s);
2073 /* double solver_params[] = {Lscale, beta, kappa, int_const, 0.0};
2074  double a0 = solve(a1,
2075  0.0,
2076  1e-9,
2077  &locked_unsat_eq,
2078  &locked_unsat_deriv,
2079  &locked_unsat_eq_deriv,
2080  static_cast<void*>(solver_params));
2081  solver_params[4] = TSTART;
2082  double astart = solve(a0,
2083  0.0,
2084  1e-9,
2085  &locked_unsat_eq,
2086  &locked_unsat_deriv,
2087  &locked_unsat_eq_deriv,
2088  static_cast<void*>(solver_params));
2089  solver_params[2] = kappa_s;
2090  solver_params[3] = int_const_s;
2091  solver_params[4] = 0;
2092  double a0_s = solve(a1,
2093  0.0,
2094  1e-9,
2095  &locked_sat_eq,
2096  &locked_sat_deriv,
2097  &locked_sat_eq_deriv,
2098  static_cast<void*>(solver_params));
2099  solver_params[4] = TSTART;
2100  double astart_s = solve(a1,
2101  0.0,
2102  1e-9,
2103  &locked_sat_eq,
2104  &locked_sat_deriv,
2105  &locked_sat_eq_deriv,
2106  static_cast<void*>(solver_params));
2107  solver_params[4] = 1.0;
2108 
2109  double w0=std::sqrt(
2110  AstroConst::G*
2111  (AstroConst::solar_mass+AstroConst::jupiter_mass)/
2112  std::pow(a0*AstroConst::solar_radius, 3))*
2113  AstroConst::day,
2114  w0_s=std::sqrt(
2115  AstroConst::G*
2116  (AstroConst::solar_mass+AstroConst::jupiter_mass)/
2117  std::pow(a0_s*AstroConst::solar_radius, 3))*
2118  AstroConst::day,
2119  w1=std::sqrt(
2120  AstroConst::G*
2121  (AstroConst::solar_mass+AstroConst::jupiter_mass)/
2122  std::pow(a1_check*AstroConst::solar_radius, 3))*
2123  AstroConst::day,
2124  w1_s=std::sqrt(
2125  AstroConst::G*
2126  (AstroConst::solar_mass+AstroConst::jupiter_mass)/
2127  std::pow(a1_check_s*AstroConst::solar_radius, 3))*
2128  AstroConst::day;*/
2129  std::valarray<double> a_transform_coef(0.0, 6),
2130  t_coef(2),
2131  t_coef_s(2);
2132  a_transform_coef[5] = Lscale / 10.0;
2133  a_transform_coef[3] = -beta / 2.0;
2134  t_coef[0] = int_const;
2135  t_coef[1] = -kappa;
2136  t_coef_s[0] = int_const_s;
2137  t_coef_s[1] = -kappa_s;
2138  std::valarray<double> Lconv_term1_coef(0.0, 2),
2139  Lconv_term2_coef(0.0, 3),
2140  identity_coef(0.0, 2),
2141  a_term1_coef_s(0.0, 3),
2142  Lconv_term1_coef_s(0.0, 2),
2143  Lc_beta_coef_s(0.0, 2);
2144  Lconv_term1_coef[1] = 1.0 / beta * std::pow(10.0 / Lscale,
2145  3.0 / 10.0);
2146  Lconv_term2_coef[2] = -2.0 / std::pow(beta, 3);
2147  identity_coef[1] = 1;
2148  a_term1_coef_s[2] = Lscale / 4.0;
2149  Lconv_term1_coef_s[1] = 1.0 / beta * std::pow(4.0 / Lscale,
2150  3.0 / 4.0);
2151  Lc_beta_coef_s[1] = 1.0 / beta;
2152 
2153  StellarEvolution::PolynomialEvolutionQuantity
2154  a_transform(a_transform_coef, 0.0, Core::Inf),
2155  a_transform1_s(a_term1_coef_s, 0.0, Core::Inf),
2156  transformed_a_evol(t_coef, TSTART, 1.0),
2157  transformed_a_evol_s(t_coef_s, TSTART, 1.0),
2158  Lconv_term1_poly(Lconv_term1_coef, 0.0, Core::Inf),
2159  Lconv_term2_poly(Lconv_term2_coef, 0.0, Core::Inf),
2160  Lconv_term1_poly_s(Lconv_term1_coef_s, 0.0, Core::Inf),
2161  Lc_beta_s(Lc_beta_coef_s, 0.0, Core::Inf),
2162  identity(identity_coef, 0.0, Core::Inf);
2163  FunctionToPower L_transform1(&Lconv_term1_poly, -10.0 / 3.0),
2164  L_transform2(&Lconv_term2_poly, -1.0),
2165  L_transform1_s(&Lconv_term1_poly_s, -4.0 / 3.0);
2166  LogFunction log_a(&identity),
2167  log_Lc_beta_s(&Lc_beta_s);
2168  ScaledFunction a_transform2_s(&log_a, -1.5 * beta),
2169  L_transform2_s(&log_Lc_beta_s, beta);
2170  FuncPlusFunc L_transform(&L_transform1, &L_transform2),
2171  a_transform_s(&a_transform1_s, &a_transform2_s),
2172  L_transform_s(&L_transform1_s, &L_transform2_s);
2173 #if 0
2174  std::valarray<double> initial_orbit(2);
2175 
2176  initial_orbit[0]=astart;
2177  initial_orbit[1]=0.0;
2178  std::cout.precision(16);
2179  std::cout.setf(std::ios_base::scientific);
2180  Star star_not_saturated_wind_no_coupling(1.0, 0.0, Kwind, wsat, Inf,
2181  0.0, 0.0, 0.0, no_evol);
2182  Planet planet1(&star_not_saturated_wind_no_coupling, 1.0, 1.0, 1.0);
2183  StellarSystem system1(&star_not_saturated_wind_no_coupling, &planet1);
2184  OrbitSolver solver(tstart, 1.0, 1e-9);
2185  solver(system1, Inf, 0.0, a0, tstart, LOCKED_TO_PLANET, initial_orbit, true);
2186  TransformedSolution to_check(a_transform, L_transform, identity, -Inf);
2187  to_check(solver);
2188  test_solution(to_check, transformed_a_evol, transformed_a_evol,
2189  zero_func, tstart, 1.0, expected_mode);
2190 
2191  initial_orbit[0]=astart_s;
2192  initial_orbit[1]=0.0;
2193  Star star_saturated_wind_no_coupling(1.0, 0.0, Kwind_s, wsat_s, Inf,
2194  0.0, 0.0, 0.0, no_evol);
2195  Planet planet2(&star_saturated_wind_no_coupling, 1.0, 1.0, 1.0);
2196  StellarSystem system2(&star_saturated_wind_no_coupling, &planet2);
2197  solver(system2, Inf, 0.0, a0_s, tstart, LOCKED_TO_PLANET,
2198  initial_orbit, true);
2199  TransformedSolution to_check_s(a_transform_s, L_transform_s,
2200  identity, -Inf);
2201  to_check_s(solver);
2202  test_solution(to_check_s, transformed_a_evol_s, transformed_a_evol_s,
2203  zero_func, tstart, 1.0, expected_mode);
2204 #endif
2205  delete no_evol;
2206  } catch (Core::Error::General &ex) {
2207  TEST_ASSERT_MSG(false, (std::string("Unexpected exception thrown: ")
2208  +
2209  ex.what()+": "+ex.get_message()).c_str());
2210  } catch (std::exception &ex) {
2211  TEST_ASSERT_MSG(false, (std::string("Unexpected exception thrown: ")
2212  +
2213  ex.what()).c_str());
2214  }
2215  }
2216 
2217  void test_OrbitSolver::test_disklocked_to_locked_to_noplanet()
2218  {
2219  try {
2220  const double Ic = 0.001,
2221  Kwind = 2e-4,
2222  tdisk = 1,
2223  tfinal = 2;
2224  StellarEvolution::MockStellarEvolution *
2225  no_evol = StellarEvolution::make_no_evolution(1.0, Ic);
2226 
2227  double afinal = 1.75,
2228  Lscale = (Core::AstroConst::jupiter_mass
2229  /
2230  std::pow(Core::AstroConst::solar_radius, 1.5)
2231  *
2232  std::sqrt(
2233  Core::AstroConst::G
2234  /
2235  (
2236  Core::AstroConst::jupiter_mass
2237  +
2238  Core::AstroConst::solar_mass
2239  )
2240  )
2241  *
2242  Core::AstroConst::day),
2243  beta = (
2244  Ic
2245  *
2246  std::sqrt(
2247  Core::AstroConst::G
2248  *
2249  (
2250  Core::AstroConst::solar_mass
2251  +
2252  Core::AstroConst::jupiter_mass
2253  )
2254  )
2255  *
2256  Core::AstroConst::day
2257  /
2258  std::pow(Core::AstroConst::solar_radius, 1.5)
2259  ),
2260  wsat = 1e100,
2261  kappa = (
2262  Kwind
2263  *
2264  std::pow(
2265  Core::AstroConst::G
2266  *
2267  (
2268  Core::AstroConst::solar_mass
2269  +
2270  Core::AstroConst::jupiter_mass
2271  )
2272  /
2273  std::pow(Core::AstroConst::solar_radius, 3),
2274  1.5
2275  )
2276  *
2277  std::pow(Core::AstroConst::day, 3)
2278  ),
2279  int_const = (Lscale / 10.0 * std::pow(afinal, 5)
2280  -
2281  beta / 2.0 * std::pow(afinal, 3)
2282  +
2283  kappa * tfinal);
2284  double solver_params[]={Lscale, beta, kappa, int_const, tdisk};
2285  double ainitial = solve(2.0 * afinal,
2286  0.0,
2287  1e-9,
2288  &locked_unsat_eq,
2289  &locked_unsat_deriv,
2290  &locked_unsat_eq_deriv,
2291  static_cast<void*>(solver_params)),
2292  wdisk = (
2293  std::sqrt(
2294  Core::AstroConst::G
2295  *
2296  (
2297  Core::AstroConst::solar_mass
2298  +
2299  Core::AstroConst::jupiter_mass
2300  )
2301  /
2302  std::pow(ainitial * Core::AstroConst::solar_radius, 3)
2303  )
2304  *
2305  Core::AstroConst::day
2306  );
2307 
2308  __star = make_const_lag_star(
2309  *no_evol,
2310  Kwind,
2311  wsat,
2312  Core::Inf,//core-env coupling timescale
2313  lag_from_lgQ(0.1, (Core::AstroConst::jupiter_mass
2314  /
2315  Core::AstroConst::solar_mass))
2316  );
2317  Planet::Planet planet(Mjup_to_Msun, Rjup_to_Rsun);
2318  planet.configure(true, //init
2319  tdisk, //age
2320  __star->mass(), //mass
2321  ainitial * (1.0 - 1e-14),//planet formation semimajor
2322  0.0, //eccentricity
2323  &zero, //spin angmom
2324  NULL, //inclination
2325  NULL, //periapsis
2326  false, //locked surface
2327  true, //zero outer inclination
2328  true);//zero outer periapsis
2329 
2330  __system = new Evolve::DiskBinarySystem(
2331  *__star,
2332  planet,
2333  ainitial * (1.0 - 1e-14), //semimajor
2334  0.0, //eccentricity
2335  0, //inclination
2336  wdisk, //Wdisk
2337  tdisk, //disk dissipation age
2338  tdisk //planet formation age
2339  );
2340 
2341  __system->configure(true, //init
2342  TSTART,
2343  Core::NaN, //semimajor
2344  Core::NaN, //eccentricity
2345  &zero, //spin angmom
2346  NULL, //inclination
2347  NULL, //periapsis
2348  Core::LOCKED_SURFACE_SPIN);
2349 
2350 
2351  double adestr = __system->minimum_separation(),
2352  tdestr = ((beta / 2.0 * std::pow(adestr, 3)
2353  +
2354  int_const
2355  -
2356  Lscale / 10.0 * std::pow(adestr, 5)) / kappa),
2357  Lc_at_destr = (beta / std::pow(adestr, 1.5)
2358  +
2359  Lscale * std::sqrt(adestr));
2360 
2361 
2362  std::valarray<double> a_transform_coef(0.0, 6), t_coef(2);
2363  a_transform_coef[5] = Lscale / 10.0;
2364  a_transform_coef[3] = -beta / 2.0;
2365  t_coef[0] = int_const;
2366  t_coef[1] = -kappa;
2367  std::valarray<double> Lconv_term1_coef(0.0, 2),
2368  Lconv_term2_coef(0.0, 3),
2369  identity_coef(0.0, 2),
2370  noplanet_Lconv_m2_coef(2);
2371  Lconv_term1_coef[1] = 1.0 / beta * std::pow(10.0 / Lscale, 0.3);
2372  Lconv_term2_coef[2] = -2.0 / std::pow(beta, 3);
2373  noplanet_Lconv_m2_coef[0] = (
2374  std::pow(Lc_at_destr, -2)
2375  -
2376  2.0 * Kwind * tdestr / std::pow(Ic, 3)
2377  );
2378  noplanet_Lconv_m2_coef[1] = 2.0 * Kwind / std::pow(Ic, 3);
2379  identity_coef[1] = 1;
2380  StellarEvolution::PolynomialEvolutionQuantity
2381  identity(identity_coef, -Core::Inf, Core::Inf),
2382  disk_nan_evol(std::valarray<double>(Core::NaN, 2),
2383  TSTART,
2384  tdisk),
2385  locked_a_transform(a_transform_coef, -Core::Inf, Core::Inf),
2386  locked_aLconv_evol(t_coef, tdisk, tdestr),
2387  locked_e_evol(std::valarray<double>(), tdisk, tdestr),
2388  noplanet_nan_evol(std::valarray<double>(Core::NaN, 2),
2389  tdestr,
2390  tfinal),
2391  disk_Lconv_evol(std::valarray<double>(wdisk * Ic, 1),
2392  TSTART,
2393  tdisk),
2394  noplanet_Lconv_m2_evol(noplanet_Lconv_m2_coef, tdestr, tfinal),
2395  Lconv_term1_poly(Lconv_term1_coef, -Core::Inf, Core::Inf),
2396  Lconv_term2_poly(Lconv_term2_coef, -Core::Inf, Core::Inf),
2397  Lrad_evol(std::valarray<double>(), TSTART, tfinal);
2398  FunctionToPower L_transform1(&Lconv_term1_poly, -10.0 / 3.0),
2399  L_transform2(&Lconv_term2_poly, -1.0),
2400  noplanet_Lconv_evol(&noplanet_Lconv_m2_evol, -0.5);
2401  LogFunction log_a(&identity);
2402  FuncPlusFunc locked_Lconv_transform(&L_transform1, &L_transform2);
2403  PiecewiseFunction a_evol, e_evol, Lconv_evol;
2404  a_evol.add_piece(&disk_nan_evol);
2405  a_evol.add_piece(&locked_aLconv_evol);
2406  a_evol.add_piece(&noplanet_nan_evol);
2407  e_evol.add_piece(&disk_nan_evol);
2408  e_evol.add_piece(&locked_e_evol);
2409  e_evol.add_piece(&noplanet_nan_evol);
2410  Lconv_evol.add_piece(&disk_Lconv_evol);
2411  Lconv_evol.add_piece(&locked_aLconv_evol);
2412  Lconv_evol.add_piece(&noplanet_Lconv_evol);
2413 
2414  std::vector< const Core::OneArgumentDiffFunction * >
2415  transformations(NUM_REAL_QUANTITIES - 1, &identity);
2416  TransformedSolution to_check(transformations, TSTART);
2417  transformations[SEMIMAJOR] = &locked_a_transform;
2418  transformations[CONV_ANGMOM] = &locked_Lconv_transform;
2419  to_check.add_transformation(transformations, tdisk);
2420  transformations[SEMIMAJOR] = &identity;
2421  transformations[CONV_ANGMOM] = &identity;
2422  to_check.add_transformation(transformations, tdestr);
2423 
2424  ExpectedEvolutionMode<Core::EvolModeType> expected_evol_mode;
2425  expected_evol_mode.add_break(TSTART, Core::LOCKED_SURFACE_SPIN);
2426  expected_evol_mode.add_break(tdisk, Core::BINARY);
2427  expected_evol_mode.add_break(tdestr, Core::SINGLE);
2428 
2429  ExpectedEvolutionMode<bool> expected_wind_mode;
2430  expected_wind_mode.add_break(TSTART, false);
2431 
2432  __star->detect_saturation();
2433  __solver = new Evolve::OrbitSolver(tfinal, 1e-8);
2434  (*__solver)(*__system,
2435  (tfinal - __system->age()) / 10000.0, //time step
2436  std::list<double>()); //no required ages*/
2437 
2438  std::vector< const std::list<double> * >
2439  tabulated_evolution = get_evolution();
2440  const std::vector< const std::list<double> * > &
2441  transformed_evolution = to_check(tabulated_evolution);
2442 
2443  std::vector<const Core::OneArgumentDiffFunction *>
2444  expected_real_quantities(NUM_REAL_QUANTITIES - 1);
2445  expected_real_quantities[SEMIMAJOR] = &a_evol;
2446  expected_real_quantities[ECCENTRICITY] = &e_evol;
2447  expected_real_quantities[CONV_INCLINATION] = &zero_func;
2448  expected_real_quantities[RAD_INCLINATION] = &zero_func;
2449  expected_real_quantities[CONV_PERIAPSIS] = &zero_func;
2450  expected_real_quantities[RAD_PERIAPSIS] = &zero_func;
2451  expected_real_quantities[CONV_ANGMOM] = &Lconv_evol;
2452  expected_real_quantities[RAD_ANGMOM] = &Lrad_evol;
2453 
2454  test_solution(transformed_evolution,
2455  expected_real_quantities,
2456  expected_evol_mode,
2457  expected_wind_mode,
2458  TSTART,
2459  tfinal);
2460  delete no_evol;
2461 
2462 #if 0
2463  Star star_not_saturated_wind_no_coupling(
2464  1.0,//M*
2465  1.0e-10,//Q*
2466  Kwind,
2467  wsat,
2468  Inf,//core-env coupling timescale
2469  0.0,//Q transition width
2470  wdisk,
2471  tdisk,
2472  no_evol);
2473  Planet planet1(&star_not_saturated_wind_no_coupling, 1.0, 1.0, 1.0);
2474  StellarSystem system1(&star_not_saturated_wind_no_coupling, &planet1);
2475  /* std::cout << std::endl << "alpha=" << Lscale
2476  << std::endl << "beta=" << beta
2477  << std::endl << "kappa=" << kappa
2478  << std::endl << "wdisk=" << wdisk << std::endl;
2479  std::cout << std::endl << "ainitial=" << ainitial << std::endl;
2480  std::cout << std::endl << "Destruction a=" << adestr
2481  << ", t=" << tdestr << ", " << ", Lc=" << Lc_at_destr
2482  << ", no planet time="
2483  << (Lscale*(std::pow(adestr, 5)-std::pow(afinal, 5))/10.0 -
2484  beta*(std::pow(adestr, 3)-std::pow(afinal, 3))/2.0)/kappa
2485  << std::endl;*/
2486 
2487  OrbitSolver solver(tstart, tfinal, 1e-8);
2488  solver(system1,
2489  Inf,//Max step
2490  0.0,//Planet formation age
2491  ainitial/AU_Rsun,//planet formation semimajor
2492  tstart);//Start age
2493 
2494  test_solution(to_check,
2495  a_evol,
2496  Lconv_evol,
2497  Lrad_evol,
2498  tstart,
2499  tfinal,
2500  expected_mode);
2501 #endif
2502  } catch (Core::Error::General &ex) {
2503  TEST_ASSERT_MSG(false, (std::string("Unexpected exception thrown: ")
2504  +
2505  ex.what()+": "+ex.get_message()).c_str());
2506  } catch (std::exception &ex) {
2507  TEST_ASSERT_MSG(false, (std::string("Unexpected exception thrown: ")
2508  +
2509  ex.what()).c_str());
2510  }
2511 
2512  }
2513 
2514  void test_OrbitSolver::test_disklocked_to_fast_to_noplanet()
2515  {
2516  try {
2517  const double TDISK = 1,
2518  TDESTR = 2,
2519  WDISK = 0.0,
2520  I_CONV = 1,
2521  L_SCALE = (
2522  -Core::AstroConst::jupiter_mass
2523  /
2524  std::pow(Core::AstroConst::solar_radius, 1.5)
2525  *
2526  std::sqrt(
2527  Core::AstroConst::G
2528  /
2529  (
2530  Core::AstroConst::jupiter_mass
2531  +
2532  Core::AstroConst::solar_mass
2533  )
2534  )
2535  *
2536  Core::AstroConst::day
2537  );
2538 
2539  StellarEvolution::MockStellarEvolution *
2540  no_evol = StellarEvolution::make_no_evolution(1.0, I_CONV);
2541 
2542  double lgQ = 8,
2543  alpha = (
2544  -4.5
2545  *
2546  std::sqrt(
2547  Core::AstroConst::G
2548  /
2549  (
2550  Core::AstroConst::solar_radius
2551  *
2552  Core::AstroConst::solar_mass
2553  )
2554  )
2555  *
2556  Core::AstroConst::jupiter_mass / std::pow(10.0, lgQ)
2557  *
2558  Core::AstroConst::Gyr / Core::AstroConst::solar_radius
2559  );
2560 
2561  __star = make_const_lag_star(
2562  *no_evol,
2563  0.0,
2564  1.0,
2565  Core::Inf,
2566  lag_from_lgQ(lgQ,
2567  (Core::AstroConst::jupiter_mass
2568  /
2569  Core::AstroConst::solar_mass))
2570  );
2571 
2572  double adestr = (
2573  2.44
2574  *
2575  std::pow(
2576  (
2577  __star->mass() * Core::AstroConst::solar_mass
2578  /
2579  Core::AstroConst::jupiter_mass
2580  ),
2581  1.0 / 3.0
2582  )
2583  *
2584  Core::AstroConst::jupiter_radius
2585  /
2586  Core::AstroConst::solar_radius
2587  );
2588  double a6p5_offset = (std::pow(adestr, 6.5)
2589  -
2590  6.5 * alpha * TDESTR),
2591  a_formation = std::pow(a6p5_offset + 6.5 * alpha * TDISK,
2592  1.0 / 6.5);
2593 
2594  evolve(WDISK,
2595  TDISK,
2596  a_formation,
2597  &zero);//Initial L*
2598 
2599  ExpectedEvolutionMode<Core::EvolModeType> expected_mode;
2600  expected_mode.add_break(TSTART, Core::LOCKED_SURFACE_SPIN);
2601  expected_mode.add_break(TDISK, Core::BINARY);
2602  expected_mode.add_break(TDESTR, Core::SINGLE);
2603 
2604  ExpectedEvolutionMode<bool> unsat_wind_mode;
2605  unsat_wind_mode.add_break(TSTART, false);
2606 
2607  std::valarray<double> a6p5_poly_coef(2);
2608  a6p5_poly_coef[0] = a6p5_offset;
2609  a6p5_poly_coef[1] = 6.5 * alpha;
2610  StellarEvolution::PolynomialEvolutionQuantity
2611  a6p5_evol(a6p5_poly_coef, TDISK, TDESTR),
2612  Lconv_disk(std::valarray<double>(I_CONV*WDISK, 1),
2613  TSTART,
2614  TDISK),
2615  Lconv_noplanet(
2616  std::valarray<double>(
2617  -L_SCALE * std::sqrt(a_formation) + I_CONV * WDISK,
2618  1
2619  ),
2620  TDESTR,
2621  MAX_AGE
2622  ),
2623  nan_disk(std::valarray<double>(Core::NaN, 2), TSTART, TDISK),
2624  nan_noplanet(std::valarray<double>(Core::NaN, 2),
2625  TDESTR,
2626  MAX_AGE),
2627  e_fast(std::valarray<double>(0.0, 1),
2628  TDISK,
2629  TDESTR),
2630  Lrad_evol(std::valarray<double>(), TSTART, MAX_AGE);
2631  FunctionToPower a_fast(&a6p5_evol, 1.0 / 6.5),
2632  sqrta_evol(&a6p5_evol, 1.0 / 13.0);
2633  ExponentialPlusFunc Lconv_unscaled(
2634  &sqrta_evol,
2635  I_CONV * WDISK / L_SCALE - std::sqrt(a_formation),
2636  0
2637  );
2638  ScaledFunction Lconv_fast(&Lconv_unscaled, L_SCALE);
2639 
2640  std::vector<const Core::OneArgumentDiffFunction *>
2641  expected_real_quantities(NUM_REAL_QUANTITIES - 1);
2642 
2643  PiecewiseFunction a_evol, e_evol, conv_angmom_evol;
2644 
2645  a_evol.add_piece(&nan_disk);
2646  a_evol.add_piece(&a_fast);
2647  a_evol.add_piece(&nan_noplanet);
2648 
2649  e_evol.add_piece(&nan_disk);
2650  e_evol.add_piece(&e_fast);
2651  e_evol.add_piece(&nan_noplanet);
2652 
2653  conv_angmom_evol.add_piece(&Lconv_disk);
2654  conv_angmom_evol.add_piece(&Lconv_fast);
2655  conv_angmom_evol.add_piece(
2656  &Lconv_noplanet
2657  );
2658 
2659  expected_real_quantities[SEMIMAJOR] = &a_evol;
2660  expected_real_quantities[ECCENTRICITY] = &e_evol;
2661  expected_real_quantities[CONV_INCLINATION] = &zero_func;
2662  expected_real_quantities[RAD_INCLINATION] = &zero_func;
2663  expected_real_quantities[CONV_PERIAPSIS] = &zero_func;
2664  expected_real_quantities[RAD_PERIAPSIS] = &zero_func;
2665  expected_real_quantities[RAD_ANGMOM] = &zero_func;
2666  expected_real_quantities[CONV_ANGMOM] = &conv_angmom_evol;
2667  expected_real_quantities[RAD_ANGMOM] = &zero_func;
2668 
2669  test_solution(get_evolution(),
2670  expected_real_quantities,
2671  expected_mode,
2672  unsat_wind_mode,
2673  TSTART,
2674  MAX_AGE);
2675 
2676  delete no_evol;
2677  } catch (Core::Error::General &ex) {
2678  TEST_ASSERT_MSG(false, (std::string("Unexpected exception thrown: ")+
2679  ex.what()+": "+ex.get_message()).c_str());
2680  } catch (std::exception &ex) {
2681  TEST_ASSERT_MSG(false, (std::string("Unexpected exception thrown: ")+
2682  ex.what()).c_str());
2683  }
2684  }
2685 
2686  void test_OrbitSolver::test_disklocked_to_fast_to_locked()
2687  {
2688  try {
2689  const double lgQ = 8,
2690  tdisk = 1,
2691  async = 2.5,
2692  tsync = 2.0,
2693  tend = 3;
2694 
2695  std::vector<const Core::OneArgumentDiffFunction *>
2696  expected_real_quantities
2697  =
2698  calculate_expected_disklocked_to_fast_to_locked(
2699  lgQ,
2700  tdisk,
2701  async,
2702  tsync,
2703  tend
2704  );
2705 
2706  double wsync = Core::orbital_angular_velocity(1.0,
2707  Mjup_to_Msun,
2708  async),
2709  Lconv_sync = (*expected_real_quantities[CONV_ANGMOM])(
2710  (tsync + tend) / 2.0
2711  ),
2712  Iconv = Lconv_sync / wsync,
2713  Ldisk = (*expected_real_quantities[CONV_ANGMOM])(
2714  (TSTART + tdisk) / 2.0
2715  ),
2716  wdisk = Ldisk / Iconv,
2717  a_formation = (*expected_real_quantities[SEMIMAJOR])(tdisk),
2718  phase_lag = lag_from_lgQ(lgQ,
2719  (Core::AstroConst::jupiter_mass
2720  /
2721  Core::AstroConst::solar_mass));
2722 
2723  StellarEvolution::MockStellarEvolution *
2724  no_evol = StellarEvolution::make_no_evolution(1.0, Iconv);
2725 
2726  __star = make_const_lag_star(
2727  *no_evol,//evolution
2728  0.0,//Kwind
2729  100.0,//wsat
2730  Core::Inf,//tcoup
2731  phase_lag
2732  );
2733 
2734  ExpectedEvolutionMode<Core::EvolModeType> expected_evol_mode;
2735  expected_evol_mode.add_break(TSTART, Core::LOCKED_SURFACE_SPIN);
2736  expected_evol_mode.add_break(tdisk, Core::BINARY);
2737 
2738  ExpectedEvolutionMode<bool> expected_wind_mode;
2739  expected_wind_mode.add_break(TSTART, false);
2740 
2741  while(true) {
2742  evolve(wdisk,//wdisk,
2743  tdisk,//tdisk
2744  a_formation,//initial semimajor
2745  (__star ? &zero : &Ldisk),//initial L*
2746  0.0,//initial inclination
2747  1.0,//planet_mass
2748  Core::NaN,//time of secondary formation
2749  tend,//max evolution age
2750  0.01,//planet radius
2751  1e-7);//precision
2752 
2753  test_solution(get_evolution(),
2754  expected_real_quantities,
2755  expected_evol_mode,
2756  expected_wind_mode,
2757  (__star ? TSTART : tdisk),
2758  tend);
2759  if(__star) {
2760  delete __star;
2761  __star = NULL;
2762 
2763  __primary_planet = new Planet::Planet(1.0, 1.0, Iconv);
2764  __primary_planet->zone().setup(
2765  std::vector<double>(),//Wtide breaks
2766  std::vector<double>(),//W* breaks
2767  std::vector<double>(1, 0.0),//Wtide pow.
2768  std::vector<double>(1, 0.0),//W* pow.
2769  phase_lag
2770  );
2771  } else
2772  break;
2773  }
2774  delete no_evol;
2775 
2776  } catch (Core::Error::General &ex) {
2777  TEST_ASSERT_MSG(
2778  false,
2779  (
2780  std::string("Unexpected exception thrown: ")
2781  +
2782  ex.what() + ": " + ex.get_message()).c_str()
2783  );
2784  } catch (std::exception &ex) {
2785  TEST_ASSERT_MSG(
2786  false,
2787  (
2788  std::string("Unexpected exception thrown: ")
2789  +
2790  ex.what()
2791  ).c_str()
2792  );
2793  }
2794  }
2795 
2796  void test_OrbitSolver::test_disklocked_to_locked_to_fast()
2797  {
2798  try {
2799  const double
2800  a0=3.2,
2801  abreak=3.0,
2802  adeath=2.0,
2803  tdisk=1,
2804  tbreak=2,
2805  tdeath=3,
2806  tend=4,
2807  Rp = (
2808  adeath
2809  *
2810  Core::AstroConst::solar_radius
2811  /
2812  Core::AstroConst::jupiter_radius
2813  /
2814  (
2815  2.44 * std::pow(Core::AstroConst::solar_mass
2816  /
2817  Core::AstroConst::jupiter_mass,
2818  1.0 / 3.0)
2819  )
2820  ),
2821  beta = (
2822  std::sqrt(
2823  Core::AstroConst::G
2824  *
2825  (
2826  Core::AstroConst::solar_mass
2827  +
2828  Core::AstroConst::jupiter_mass
2829  )
2830  )
2831  *
2832  Core::AstroConst::day
2833  /
2834  std::pow(Core::AstroConst::solar_radius, 1.5)
2835  ),
2836  wdisk = beta / std::pow(a0, 1.5),
2837  gamma = (
2838  (std::pow(abreak, 6.5) - std::pow(adeath, 6.5))
2839  /
2840  (tdeath - tbreak)
2841  ),
2842  Q = (
2843  9.0 * 13.0 / 4.0
2844  *
2845  std::sqrt(Core::AstroConst::G
2846  /
2847  Core::AstroConst::solar_mass)
2848  *
2849  Core::AstroConst::jupiter_mass
2850  *
2851  Core::AstroConst::Gyr
2852  /
2853  (gamma * std::pow(Core::AstroConst::solar_radius, 1.5))
2854  ),
2855  alphaL = (
2856  Core::AstroConst::jupiter_mass
2857  /
2858  std::pow(Core::AstroConst::solar_radius, 1.5)
2859  *
2860  std::sqrt(Core::AstroConst::G
2861  /
2862  (
2863  Core::AstroConst::jupiter_mass
2864  +
2865  Core::AstroConst::solar_mass
2866  )
2867  )
2868  *
2869  Core::AstroConst::day
2870  ),
2871  Ic = (
2872  (
2873  (
2874  (std::pow(a0, 5) - std::pow(abreak, 5))
2875  *
2876  std::pow(abreak, 3.5)
2877  *
2878  6.5
2879  )
2880  -
2881  5.0 * gamma * abreak * abreak
2882  )
2883  /
2884  (
2885  (
2886  (std::pow(a0, 3) - std::pow(abreak, 3))
2887  *
2888  std::pow(abreak, 3.5) * 6.5
2889  )
2890  -
2891  3.0 * gamma
2892  )
2893  *
2894  alphaL
2895  /
2896  (5.0 * beta)
2897  ),
2898  kappa = (
2899  (std::pow(a0, 5) - std::pow(abreak, 5)) * alphaL / 10.0
2900  -
2901  (std::pow(a0, 3) - std::pow(abreak, 3)) * beta * Ic / 2.0
2902  ),
2903  Kwind = (
2904  kappa
2905  /
2906  std::pow(
2907  Core::AstroConst::G
2908  *
2909  (
2910  Core::AstroConst::solar_mass
2911  +
2912  Core::AstroConst::jupiter_mass
2913  )
2914  /
2915  std::pow(Core::AstroConst::solar_radius, 3),
2916  1.5
2917  )
2918  /
2919  std::pow(Core::AstroConst::day, 3)
2920  ),
2921  locked_a_int_const = (alphaL * std::pow(a0, 5) / 10.0
2922  -
2923  beta * Ic * std::pow(a0, 3) / 2.0
2924  +
2925  kappa * tdisk);
2926 
2927  StellarEvolution::MockStellarEvolution *
2928  no_evol = StellarEvolution::make_no_evolution(1.0, Ic);
2929 
2930  __star = make_const_lag_star(
2931  *no_evol,//evolution
2932  Kwind,//Kwind
2933  100.0,//wsat
2934  Core::Inf,//tcoup
2935  lag_from_lgQ(
2936  std::log10(Q),
2937  (Core::AstroConst::jupiter_mass
2938  /
2939  Core::AstroConst::solar_mass)
2940  )
2941  );
2942  evolve(wdisk,//wdisk,
2943  tdisk,//tdisk
2944  a0 * (1.0 + 1e-14),//initial semimajor
2945  &zero,//initial L*
2946  0.0,//initial inclination
2947  1.0,//planet_mass
2948  Core::NaN,//form the planet when disk dissipates
2949  tend,//max evolution age
2950  Rp);//planet radius
2951 
2952  ExpectedEvolutionMode<Core::EvolModeType> expected_evol_mode;
2953  expected_evol_mode.add_break(TSTART, Core::LOCKED_SURFACE_SPIN);
2954  expected_evol_mode.add_break(tdisk, Core::BINARY);
2955  expected_evol_mode.add_break(tdeath, Core::SINGLE);
2956 
2957  ExpectedEvolutionMode<bool> expected_wind_mode;
2958  expected_wind_mode.add_break(TSTART, false);
2959 
2960  std::valarray<double> a_locked_transform_coef(0.0, 6),
2961  a_locked_evol_coef(2),
2962  identity_coef(0.0, 2),
2963  a6p5_fast_evol_coef(2),
2964  Lconv_locked_term1_coef(0.0, 2),
2965  Lconv_locked_term2_coef(0.0, 3);
2966  a_locked_transform_coef[5] = alphaL / 10.0;
2967  a_locked_transform_coef[3] = -beta * Ic / 2.0;
2968  a_locked_evol_coef[0] = locked_a_int_const;
2969  a_locked_evol_coef[1] = -kappa;
2970  identity_coef[1] = 1.0;
2971  a6p5_fast_evol_coef[0] = gamma * tbreak + std::pow(abreak, 6.5);
2972  a6p5_fast_evol_coef[1] = -gamma;
2973  Lconv_locked_term1_coef[1] = (1.0
2974  /
2975  (beta * Ic)
2976  *
2977  std::pow(10.0 / alphaL, 3.0 / 10.0));
2978  Lconv_locked_term2_coef[2] = -2.0 / std::pow(beta * Ic, 3);
2979 
2980  StellarEvolution::PolynomialEvolutionQuantity
2981  identity(identity_coef, -Core::Inf, Core::Inf),
2982  a_disk_transform = identity,
2983  nan_disk_evol(std::valarray<double>(Core::NaN, 2),
2984  TSTART,
2985  tdisk),
2986  a_locked_transform(a_locked_transform_coef,
2987  -Core::Inf,
2988  Core::Inf),
2989  a_locked_evol(a_locked_evol_coef, tdisk, tbreak),
2990  a_fast_transform = identity,
2991  a6p5_fast_evol(a6p5_fast_evol_coef, tbreak, tdeath),
2992  a_noplanet_transform = identity,
2993  nan_noplanet_evol(std::valarray<double>(Core::NaN, 2),
2994  tdeath,
2995  tend),
2996 
2997  Lconv_disk_transform = identity,
2998  Lconv_disk_evol(std::valarray<double>(wdisk * Ic, 1),
2999  TSTART,
3000  tdisk),
3001  Lconv_locked_term1_poly(Lconv_locked_term1_coef,
3002  -Core::Inf,
3003  Core::Inf),
3004  Lconv_locked_term2_poly(Lconv_locked_term2_coef,
3005  -Core::Inf,
3006  Core::Inf),
3007  Lconv_locked_evol = a_locked_evol,
3008  Lconv_fast_transform(std::valarray<double>(Core::NaN, 2),
3009  -Core::Inf,
3010  Core::Inf),
3011  Lconv_fast_evol(std::valarray<double>(Core::NaN, 2),
3012  tbreak,
3013  tdeath),
3014  Lconv_noplanet_transform(std::valarray<double>(Core::NaN, 2),
3015  -Core::Inf,
3016  Core::Inf),
3017  Lconv_noplanet_evol(std::valarray<double>(Core::NaN, 2),
3018  tdeath,
3019  tend),
3020  Lrad_transform = identity,
3021  Lrad_evol(std::valarray<double>(), TSTART, tend),
3022  zero_e(std::valarray<double>(), tdisk, tdeath);
3023 
3024  FunctionToPower
3025  a_fast_evol(&a6p5_fast_evol, 1.0 / 6.5),
3026  Lconv_locked_transform1(&Lconv_locked_term1_poly,
3027  -10.0 / 3.0),
3028  Lconv_locked_transform2(&Lconv_locked_term2_poly, -1.0);
3029 
3030  FuncPlusFunc Lconv_locked_transform(&Lconv_locked_transform1,
3031  &Lconv_locked_transform2);
3032 
3033  PiecewiseFunction a_evol, e_evol, Lconv_evol;
3034 
3035  a_evol.add_piece(&nan_disk_evol);
3036  a_evol.add_piece(&a_locked_evol);
3037  a_evol.add_piece(&a_fast_evol);
3038  a_evol.add_piece(&nan_noplanet_evol);
3039 
3040  e_evol.add_piece(&nan_disk_evol);
3041  e_evol.add_piece(&zero_e);
3042  e_evol.add_piece(&nan_noplanet_evol);
3043 
3044  Lconv_evol.add_piece(&Lconv_disk_evol);
3045  Lconv_evol.add_piece(&Lconv_locked_evol);
3046  Lconv_evol.add_piece(&Lconv_fast_evol);
3047  Lconv_evol.add_piece(&Lconv_noplanet_evol);
3048 
3049  std::vector< const Core::OneArgumentDiffFunction * >
3050  transformations(NUM_REAL_QUANTITIES - 1, &identity);
3051  TransformedSolution to_check(transformations, TSTART);
3052 
3053  transformations[SEMIMAJOR] = &a_disk_transform;
3054  transformations[CONV_ANGMOM] = &Lconv_disk_transform;
3055  to_check.add_transformation(transformations, TSTART);
3056 
3057  transformations[SEMIMAJOR] = &a_locked_transform;
3058  transformations[CONV_ANGMOM] = &Lconv_locked_transform;
3059  to_check.add_transformation(transformations, tdisk);
3060 
3061  transformations[SEMIMAJOR] = &a_fast_transform;
3062  transformations[CONV_ANGMOM] = &Lconv_fast_transform;
3063  to_check.add_transformation(transformations, tbreak);
3064 
3065  transformations[SEMIMAJOR] = &a_fast_transform;
3066  transformations[CONV_ANGMOM] = &Lconv_fast_transform;
3067  to_check.add_transformation(transformations, tbreak);
3068 
3069  transformations[SEMIMAJOR] = &a_noplanet_transform;
3070  transformations[CONV_ANGMOM] = &Lconv_noplanet_transform;
3071  to_check.add_transformation(transformations, tdeath);
3072 
3073  std::vector<const Core::OneArgumentDiffFunction *>
3074  expected_real_quantities(NUM_REAL_QUANTITIES - 1);
3075 
3076  std::vector< const std::list<double> * >
3077  tabulated_evolution = get_evolution();
3078 
3079  const std::vector< const std::list<double> * > &
3080  transformed_evolution = to_check(tabulated_evolution);
3081 
3082  expected_real_quantities[SEMIMAJOR] = &a_evol;
3083  expected_real_quantities[ECCENTRICITY] = &e_evol;
3084  expected_real_quantities[CONV_INCLINATION] = &zero_func;
3085  expected_real_quantities[RAD_INCLINATION] = &zero_func;
3086  expected_real_quantities[CONV_PERIAPSIS] = &zero_func;
3087  expected_real_quantities[RAD_PERIAPSIS] = &zero_func;
3088  expected_real_quantities[CONV_ANGMOM] = &Lconv_evol;
3089  expected_real_quantities[RAD_ANGMOM] = &zero_func;
3090 
3091  test_solution(transformed_evolution,
3092  expected_real_quantities,
3093  expected_evol_mode,
3094  expected_wind_mode,
3095  TSTART,
3096  tend);
3097  delete no_evol;
3098  } catch (Core::Error::General &ex) {
3099  TEST_ASSERT_MSG(
3100  false,
3101  (
3102  std::string("Unexpected exception thrown: ")
3103  +
3104  ex.what()+": "+ex.get_message()
3105  ).c_str()
3106  );
3107  } catch (std::exception &ex) {
3108  TEST_ASSERT_MSG(
3109  false,
3110  (
3111  std::string("Unexpected exception thrown: ")
3112  +
3113  ex.what()
3114  ).c_str()
3115  );
3116  }
3117 
3118  }
3119 
3120  void test_OrbitSolver::test_polar_1_0_evolution()
3121  {
3122  StellarEvolution::MockStellarEvolution *
3123  no_evol = StellarEvolution::make_no_evolution();
3124 
3125  const double TDISK = 0.1,
3126  WSTAR = 0.01,
3127  WORB = 0.1;
3128  double initial_L = 1.0;
3129  std::vector<const Core::OneArgumentDiffFunction *>
3130  expected_real_quantities = calculate_expected_polar_1_0(
3131  TDISK,
3132  WSTAR,
3133  WORB
3134  );
3135  ExpectedEvolutionMode<Core::EvolModeType> expected_evol_mode;
3136  expected_evol_mode.add_break(TSTART, Core::LOCKED_SURFACE_SPIN);
3137  expected_evol_mode.add_break(TDISK, Core::BINARY);
3138 
3139  ExpectedEvolutionMode<bool> expected_wind_mode;
3140  expected_wind_mode.add_break(TSTART, false);
3141 
3142  double semimajor = (*expected_real_quantities[SEMIMAJOR])(
3143  (TDISK + MAX_AGE) / 2.0
3144  );
3145 
3146  make_single_component_star(*no_evol,
3147  0.0,//Kw
3148  1.0,//Wsat
3149  Core::Inf,
3150  0.9 * WSTAR,
3151  1.1 * WSTAR,
3152  0.1 * WSTAR,
3153  1.0);//tcoup
3154 
3155  while(true) {
3156  evolve(WSTAR,//wdisk
3157  TDISK,//tdisk
3158  semimajor,//initial semimajor
3159  &initial_L,//initial L*
3160  M_PI / 2.0);//initial inclination
3161 
3162  test_solution(get_evolution(),
3163  expected_real_quantities,
3164  expected_evol_mode,
3165  expected_wind_mode,
3166  (__star ? TSTART : TDISK),
3167  MAX_AGE);
3168 
3169  if(__star) {
3170  delete __star;
3171  delete __system;
3172  delete __solver;
3173  __star = NULL;
3174  __system = NULL;
3175  __solver = NULL;
3176 
3177  __primary_planet = new Planet::Planet(1.0, 1.0, 1.0);
3178  set_single_component_dissipation(
3179  0.9 * WSTAR,
3180  1.1 * WSTAR,
3181  0.1 * WSTAR,
3182  1.0
3183  );
3184  initial_L = WSTAR;
3185  } else
3186  break;
3187  }
3188 
3189  delete no_evol;
3190  }
3191 
3192  void test_OrbitSolver::test_polar_2_0_evolution()
3193  {
3194  StellarEvolution::MockStellarEvolution *
3195  no_evol = StellarEvolution::make_no_evolution();
3196 
3197  const double TDISK = 0.1,
3198  WSTAR = 0.01,
3199  WORB = 0.1,
3200  PHASE_LAG = 1e-2;
3201 
3202  double lconv_decay_rate, semimajor, initial_L = 1.0;
3203 
3204  std::vector<const Core::OneArgumentDiffFunction *>
3205  expected_real_quantities = calculate_expected_polar_2_0(TDISK,
3206  WSTAR,
3207  WORB,
3208  PHASE_LAG,
3209  lconv_decay_rate,
3210  semimajor);
3211 
3212 
3213  std::valarray<double> identity_coef(0.0, 2);
3214  identity_coef[1] = 1.0;
3215  StellarEvolution::PolynomialEvolutionQuantity
3216  identity(identity_coef, -Core::Inf, Core::Inf),
3217  one_func(std::valarray<double>(1.0, 1), -Core::Inf, Core::Inf);
3218 
3219  CosFunction cos_transform(&identity);
3220  FuncPlusFunc one_plus_cos_transform(&one_func,
3221  &cos_transform);
3222 
3223 
3224  make_single_component_star(
3225  *no_evol,
3226  0.0,//Kw
3227  1.0,//Wsat
3228  Core::Inf,//tcoup
3229  2.0 * (WSTAR - lconv_decay_rate * MAX_AGE),
3230  2.0 * (WSTAR + lconv_decay_rate * MAX_AGE),
3231  0.1 * WSTAR,
3232  PHASE_LAG
3233  );
3234 
3235  while(true) {
3236  evolve(WSTAR,//wdisk
3237  TDISK,//tdisk
3238  semimajor,//initial semimajor
3239  &initial_L,//initial L*
3240  M_PI / 2.0);//initial inclination
3241 
3242  std::vector< const Core::OneArgumentDiffFunction * >
3243  transformations(NUM_REAL_QUANTITIES - 1, &identity);
3244  transformations[CONV_INCLINATION] = &one_plus_cos_transform;
3245  transformations[RAD_INCLINATION] = &one_plus_cos_transform;
3246  TransformedSolution to_check(transformations, TSTART);
3247 
3248  std::vector< const std::list<double> * >
3249  tabulated_evolution = get_evolution();
3250  const std::vector< const std::list<double> * > &
3251  transformed_evolution = to_check(tabulated_evolution);
3252 
3253  ExpectedEvolutionMode<Core::EvolModeType> expected_evol_mode;
3254  expected_evol_mode.add_break(TSTART, Core::LOCKED_SURFACE_SPIN);
3255  expected_evol_mode.add_break(TDISK, Core::BINARY);
3256 
3257  ExpectedEvolutionMode<bool> expected_wind_mode;
3258  expected_wind_mode.add_break(TSTART, false);
3259 
3260  test_solution(transformed_evolution,
3261  expected_real_quantities,
3262  expected_evol_mode,
3263  expected_wind_mode,
3264  (__star ? TSTART : TDISK),
3265  MAX_AGE);
3266 
3267  if(__star) {
3268  delete __star;
3269  delete __system;
3270  delete __solver;
3271  __star = NULL;
3272  __system = NULL;
3273  __solver = NULL;
3274 
3275  __primary_planet = new Planet::Planet(1.0, 1.0, 1.0);
3276  set_single_component_dissipation(
3277  2.0 * (WSTAR - lconv_decay_rate * MAX_AGE),
3278  2.0 * (WSTAR + lconv_decay_rate * MAX_AGE),
3279  0.1 * WSTAR,
3280  PHASE_LAG
3281  );
3282  initial_L = WSTAR;
3283  } else
3284  break;
3285  }
3286 
3287  delete no_evol;
3288  }
3289  void test_OrbitSolver::test_oblique_1_0_evolution()
3290  {
3291  StellarEvolution::MockStellarEvolution *
3292  no_evol = StellarEvolution::make_no_evolution();
3293 
3294  const double PHASE_LAG = 0.1,
3295  TDISK = 0.1,
3296  WORB = 0.1,
3297  INITIAL_INC = 0.25 * M_PI,
3298  INITIAL_WSTAR = 0.01;
3299 
3300  double initial_angmom = 1.0,
3301  min_wstar;
3302 
3303  std::vector<const Core::OneArgumentDiffFunction *>
3304  expected_real_quantities = calculate_expected_oblique_m_0(
3305  1,
3306  TDISK,
3307  WORB,
3308  INITIAL_INC,
3309  INITIAL_WSTAR,
3310  PHASE_LAG,
3311  min_wstar
3312  );
3313 
3314  double wstar_tolerance = 0.01 * (INITIAL_WSTAR - min_wstar),
3315  semimajor = (*expected_real_quantities[SEMIMAJOR])(
3316  (TDISK + MAX_AGE) / 2.0
3317  );
3318 
3319  ExpectedEvolutionMode<Core::EvolModeType> expected_evol_mode;
3320  expected_evol_mode.add_break(TSTART, Core::LOCKED_SURFACE_SPIN);
3321  expected_evol_mode.add_break(TDISK, Core::BINARY);
3322 
3323  ExpectedEvolutionMode<bool> expected_wind_mode;
3324  expected_wind_mode.add_break(TSTART, false);
3325 
3326  make_single_component_star(
3327  *no_evol,
3328  0.0,//Kw
3329  1.0,//Wsat
3330  Core::Inf,//tcoup
3331  min_wstar - wstar_tolerance,
3332  INITIAL_WSTAR + wstar_tolerance,
3333  0.1 * INITIAL_WSTAR,
3334  PHASE_LAG
3335  );
3336 
3337  while(true) {
3338  evolve(INITIAL_WSTAR,//wdisk
3339  TDISK,//tdisk
3340  semimajor,//initial semimajor
3341  &initial_angmom,//initial L*
3342  INITIAL_INC);//initial inclination
3343 
3344  test_solution(get_evolution(),
3345  expected_real_quantities,
3346  expected_evol_mode,
3347  expected_wind_mode,
3348  (__star ? TSTART : TDISK),
3349  MAX_AGE);
3350 
3351  if(__star) {
3352  delete __star;
3353  delete __system;
3354  delete __solver;
3355  __star = NULL;
3356  __system = NULL;
3357  __solver = NULL;
3358 
3359  __primary_planet = new Planet::Planet(1.0, 1.0, 1.0);
3360  set_single_component_dissipation(
3361  min_wstar - wstar_tolerance,
3362  INITIAL_WSTAR + wstar_tolerance,
3363  0.1 * INITIAL_WSTAR,
3364  PHASE_LAG
3365  );
3366 
3367  initial_angmom = INITIAL_WSTAR;
3368  } else
3369  break;
3370  }
3371  }
3372 
3373  void test_OrbitSolver::test_oblique_2_0_evolution()
3374  {
3375  StellarEvolution::MockStellarEvolution *
3376  no_evol = StellarEvolution::make_no_evolution();
3377 
3378  const double PHASE_LAG = 0.1,
3379  TDISK = 0.1,
3380  WORB = 0.1,
3381  INITIAL_INC = 0.25 * M_PI,
3382  INITIAL_WSTAR = 0.01;
3383 
3384  double initial_angmom = 1.0,
3385  min_wstar;
3386 
3387  std::vector<const Core::OneArgumentDiffFunction *>
3388  expected_real_quantities = calculate_expected_oblique_m_0(
3389  2,
3390  TDISK,
3391  WORB,
3392  INITIAL_INC,
3393  INITIAL_WSTAR,
3394  PHASE_LAG,
3395  min_wstar
3396  );
3397 
3398  double wstar_tolerance = 0.01 * (INITIAL_WSTAR - min_wstar),
3399  semimajor = (*expected_real_quantities[SEMIMAJOR])(
3400  (TDISK + MAX_AGE) / 2.0
3401  );
3402 
3403  ExpectedEvolutionMode<Core::EvolModeType> expected_evol_mode;
3404  expected_evol_mode.add_break(TSTART, Core::LOCKED_SURFACE_SPIN);
3405  expected_evol_mode.add_break(TDISK, Core::BINARY);
3406 
3407  ExpectedEvolutionMode<bool> expected_wind_mode;
3408  expected_wind_mode.add_break(TSTART, false);
3409 
3410  make_single_component_star(
3411  *no_evol,
3412  0.0,//Kw
3413  1.0,//Wsat
3414  Core::Inf,//tcoup
3415  2.0 * min_wstar - wstar_tolerance,
3416  2.0 * INITIAL_WSTAR + wstar_tolerance,
3417  0.1 * INITIAL_WSTAR,
3418  PHASE_LAG
3419  );
3420 
3421  while(true) {
3422  evolve(INITIAL_WSTAR,//wdisk
3423  TDISK,//tdisk
3424  semimajor,//initial semimajor
3425  &initial_angmom,//initial L*
3426  INITIAL_INC);//initial inclination
3427 
3428  test_solution(get_evolution(),
3429  expected_real_quantities,
3430  expected_evol_mode,
3431  expected_wind_mode,
3432  (__star ? TSTART : TDISK),
3433  MAX_AGE);
3434  if(__star) {
3435  delete __star;
3436  delete __system;
3437  delete __solver;
3438  __star = NULL;
3439  __system = NULL;
3440  __solver = NULL;
3441 
3442  __primary_planet = new Planet::Planet(1.0, 1.0, 1.0);
3443  set_single_component_dissipation(
3444  2.0 * min_wstar - wstar_tolerance,
3445  2.0 * INITIAL_WSTAR + wstar_tolerance,
3446  0.1 * INITIAL_WSTAR,
3447  PHASE_LAG
3448  );
3449 
3450  initial_angmom = INITIAL_WSTAR;
3451  } else
3452  break;
3453  }
3454 
3455  delete no_evol;
3456 
3457  }
3458 
3459  void test_OrbitSolver::setup()
3460  {
3461  assert(__star == NULL);
3462  assert(__primary_planet == NULL);
3463  assert(__system == NULL);
3464  assert(__solver == NULL);
3465  assert(__temp_functions.empty());
3466  }
3467 
3468  void test_OrbitSolver::tear_down()
3469  {
3470  if(__star) {
3471  std::cout << "Deleting star" << std::endl;
3472  delete __star;
3473  __star = NULL;
3474  }
3475  if(__primary_planet) {
3476  std::cout << "Deleting primary planet" << std::endl;
3477  delete __primary_planet;
3478  __primary_planet = NULL;
3479  }
3480  if(__system) {
3481  std::cout << "Deleting system" << std::endl;
3482  delete __system;
3483  __system = NULL;
3484  }
3485  if(__solver) {
3486  std::cout << "Deleting solver" << std::endl;
3487  delete __solver;
3488  __solver = NULL;
3489  }
3490  for(
3491  std::vector< const Core::OneArgumentDiffFunction* >::iterator
3492  temp_func_i = __temp_functions.begin();
3493  temp_func_i != __temp_functions.end();
3494  ++temp_func_i
3495  ) {
3496  std::cout << "Deleting function" << std::endl;
3497  delete *temp_func_i;
3498  }
3499  __temp_functions.clear();
3500  }
3501 
3502  test_OrbitSolver::test_OrbitSolver() :
3503  __solver(NULL),
3504  __system(NULL),
3505  __star(NULL),
3506  __primary_planet(NULL)
3507  {
3508  TEST_ADD(test_OrbitSolver::test_disk_locked_no_stellar_evolution);
3509  TEST_ADD(test_OrbitSolver::test_disk_locked_with_stellar_evolution);
3510  TEST_ADD(test_OrbitSolver::test_no_planet_evolution);
3511  TEST_ADD(test_OrbitSolver::test_unlocked_evolution);
3512 // TEST_ADD(test_OrbitSolver::test_locked_evolution);//NOT REVIVED!!!
3513  TEST_ADD(test_OrbitSolver::test_disklocked_to_locked_to_noplanet);
3514  TEST_ADD(test_OrbitSolver::test_disklocked_to_fast_to_noplanet);
3515  TEST_ADD(test_OrbitSolver::test_disklocked_to_fast_to_locked);
3516  TEST_ADD(test_OrbitSolver::test_disklocked_to_locked_to_fast);
3517  TEST_ADD(test_OrbitSolver::test_polar_1_0_evolution);
3518  TEST_ADD(test_OrbitSolver::test_polar_2_0_evolution);
3519  TEST_ADD(test_OrbitSolver::test_oblique_1_0_evolution);
3520  TEST_ADD(test_OrbitSolver::test_oblique_2_0_evolution);
3521  }
3522 
3523 }//End Evolve namespace.
StellarEvolution::MockStellarEvolution
Implements a StellarEvolution based on polynomial evolution quantities.
Definition: PolynomialEvolution.h:126
Evolve::test_OrbitSolver::__star
Star::InterpolatedEvolutionStar * __star
The star used in the current test (NULL if primary is planet).
Definition: testOrbitSolver.h:45
MIN_AGE
const double MIN_AGE
Most tests start at this age in Gyr.
Definition: Common.h:54
InCol::TDISK
Age in Myr when disk dissipates.
Definition: IOColumns.h:62
Evolve::BinarySystem::eccentricity_evolution
const std::list< double > & eccentricity_evolution() const
The tabulated evolution of the eccentricity so far.
Definition: BinarySystem.h:1052
Evolve::locked_sat_eq
double locked_sat_eq(double a, void *params)
The equation that should be solved in order to get the semimamjor axis at a given time for the locked...
Definition: testOrbitSolver.cpp:1888
Evolve::INCLINATION
Inclination of the zone.
Definition: DissipatingZone.h:36
testOrbitSolver.h
Declares the test suite that exercises the OrbitSolver class and some other clasess necessary to acco...
OutCol::RAD_PERIAPSIS
The orbital periapsis in the reference frame of the stellar radiative zone in radians.
Definition: IOColumns.h:186
Star::SaturatingSkumanichWindBody::detect_saturation
void detect_saturation()
Sets the saturation based on the currently configured spin frequency.
Definition: SaturatingSkumanichWindBody.h:65
Core::Error::General
The base class of all exceptions.
Definition: Error.h:27
Evolve::OrbitSolver::mode_evolution
const std::list< Core::EvolModeType > & mode_evolution() const
The tabulated evolution modes so far.
Definition: OrbitSolver.h:498
OutCol::WORB
The orbital frequency in rad/day.
Definition: IOColumns.h:168
ScaledFunction
A function scaled by some constant.
Definition: PolynomialEvolution.h:594
Evolve::operator<<
std::ostream & operator<<(std::ostream &os, const ZoneEvolutionQuantities &evol_var)
More civilized output for EvolVarType variables.
Definition: DissipatingZone.cpp:7
Evolve::DissipatingBody
A base class for any body contributing to tidal dissipation.
Definition: DissipatingBody.h:39
Evolve::locked_sat_eq_deriv
void locked_sat_eq_deriv(double a, void *params, double *f, double *df)
The equation and its derivative that should be solved in order to get the semimamjor axis at a given ...
Definition: testOrbitSolver.cpp:1965
Evolve::BinarySystem::minimum_separation
virtual double minimum_separation(bool deriv=false) const
Smallest semimajor axis at which the secondary can survive for the latest system configuration.
Definition: BinarySystem.cpp:2248
Core::OneArgumentDiffFunction
A class representing a once differentiable function of a single argument.
Definition: Functions.h:104
SEMIMAJOR
SEMIMAJOR
The derivative w.r.t. the semimajor axis in AU.
Definition: DissipationQuantities.h:99
Planet::Planet
Single zone non-evolving planets with huge dissipation, so they always remain locked to the disk...
Definition: Planet.h:21
Evolve::test_OrbitSolver::test_no_planet_scenario
void test_no_planet_scenario(const StellarEvolution::Interpolator &stellar_evol, double *initial_Lstar, double windK, double wind_sat_freq, double core_env_coupling_time, std::vector< const Core::OneArgumentDiffFunction *> &expected_evolution, const ExpectedEvolutionMode< bool > &expected_wind_mode, double max_age=MAX_AGE, bool debug_mode=false)
Test a planet-less scenario computed in 3 different ways: 1) withou a planet 2) without dissipation 3...
Definition: testOrbitSolver.cpp:517
OutCol::RAD_INCLINATION
The angle between the stellar core spin and orbital angular momentum in radians.
Definition: IOColumns.h:178
Planet::Planet::zone
const Evolve::DissipatingZone & zone(unsigned zone_index) const
Returns the only zone.
Definition: Planet.h:37
OrbitSolver
struct LIB_PUBLIC OrbitSolver
Opaque struct to cast to/from Evolve::OrbitSolver.
Definition: CInterface.h:38
OutCol::CONV_INCLINATION
The angle between the stellar surface spin and orbital angular momentum in radians.
Definition: IOColumns.h:174
Core::AstroConst::jupiter_mass
const double jupiter_mass
Mass of Jupiter [kg].
Definition: AstronomicalConstants.h:27
Evolve::test_OrbitSolver::make_single_component_star
void make_single_component_star(const StellarEvolution::Interpolator &evolution, double wind_strength, double wind_sat_freq, double coupling_timescale, double min_frequency, double max_frequency, double decay_scale, double phase_lag=1.0e-5)
Create __star with constant dissipation in a range, quickly decaying outside of that.
Definition: testOrbitSolver.cpp:132
Evolve::test_OrbitSolver::__system
Evolve::DiskBinarySystem * __system
The system being evolved by the current test.
Definition: testOrbitSolver.h:42
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::BinarySystem::semimajor_evolution
const std::list< double > & semimajor_evolution() const
The tabulated evolution of the semimajor axis so far.
Definition: BinarySystem.h:1044
Star::SaturatingSkumanichWindBody::wind_saturation_evolution
const std::list< bool > & wind_saturation_evolution() const
The tabulated wind saturation states so far.
Definition: SaturatingSkumanichWindBody.h:129
Star::InterpolatedEvolutionStar::core
const EvolvingStellarCore & core() const
The core of the star.
Definition: EvolvingStar.h:100
Evolve::test_OrbitSolver::__primary_planet
Planet::Planet * __primary_planet
The primary planet in the current test (NULL if primary is star).
Definition: testOrbitSolver.h:48
Evolve::RealEvolutionQuantity
RealEvolutionQuantity
Define identifiers for the quantities whose evolution we check.
Definition: RealEvolutionQuantity.h:15
DiskBinarySystem
struct LIB_PUBLIC DiskBinarySystem
Opaque struct to cast to/from Evolve::DiskBinarySystem.
Definition: CInterface.h:35
StellarEvolution::Interpolator
A class that interpolates among stellar evolution tracks.
Definition: Interpolator.h:42
FuncPlusFunc
Represents the sum of two functions and the derivative.
Definition: PolynomialEvolution.h:438
Core::AstroConst::jupiter_radius
const double jupiter_radius
Radius of Jupiter [m].
Definition: AstronomicalConstants.h:28
Evolve::DissipatingBody::mass
double mass() const
The mass of the body (constant with age).
Definition: DissipatingBody.h:534
Evolve::DissipatingBody::configure
virtual void configure(bool initialize, double age, double companion_mass, double semimajor, double eccentricity, const double *spin_angmom, const double *inclination=NULL, const double *periapsis=NULL, bool locked_surface=false, bool zero_outer_inclination=false, bool zero_outer_periapsis=false)
Defines the orbit this body is in.
Definition: DissipatingBody.cpp:603
Evolve::test_OrbitSolver::evolve
void evolve(double wdisk, double tdisk, double initial_a, const double *initial_Lstar, double initial_incl=0.0, double secondary_mass=1.0, double tsecondary=Core::NaN, double max_age=MAX_AGE, double secondary_radius=1.0, double precision=1e-6, double max_step_factor=1e-3)
Add a planet to the given star and evolve, returning the solver.
Definition: testOrbitSolver.cpp:156
lag_from_lgQ
double lag_from_lgQ(double lgQ)
Converts lg(Q) to a tidal phase lag.
Definition: CInterface.cpp:17
Evolve::PERIAPSIS
Periapsis of the zone.
Definition: DissipatingZone.h:42
Core::AstroConst::solar_radius
const double solar_radius
Radius of the sun [m].
Definition: AstronomicalConstants.h:22
Evolve::DissipatingZone
A layer of a system body for which the tidal bulge is not exactly in phase with the tidal potential...
Definition: DissipatingZone.h:132
CosFunction
The cosine of a function.
Definition: PolynomialEvolution.h:668
StellarEvolution::PolynomialEvolutionQuantity
An EvolvingStellar quantity that uses a polynomial instead of interpolating.
Definition: PolynomialEvolution.h:28
LogFunction
The natural logarithm of a function.
Definition: PolynomialEvolution.h:632
InCol::WDISK
Stellar surface spin while disk is present in rad/day.
Definition: IOColumns.h:57
Evolve::test_OrbitSolver::set_single_component_dissipation
void set_single_component_dissipation(double min_frequency, double max_frequency, double decay_scale, double phase_lag=1.0e-5)
Set the dissipation of the primary to only a single tidal.
Definition: testOrbitSolver.cpp:84
Core::Error::General::what
virtual const char * what() const
The type of error.
Definition: Error.h:39
Core::AstroConst::Gyr
const double Gyr
Gyr [s].
Definition: AstronomicalConstants.h:24
check_diff
bool check_diff(double x, double y, double frac_tolerance, double abs_tolerance)
Returns true iff .
Definition: Common.cpp:3
Evolve::locked_unsat_deriv
double locked_unsat_deriv(double a, void *params)
The derivative of the equation that should be solved in order to get the semimamjor axis at a given t...
Definition: testOrbitSolver.cpp:1912
ECCENTRICITY
ECCENTRICITY
The derivative w.r.t. the eccentricity.
Definition: DissipationQuantities.h:88
Core::AstroConst::day
const double day
day [s]
Definition: AstronomicalConstants.h:25
Evolve::test_OrbitSolver::calculate_expected_unlocked_evolution
std::vector< const Core::OneArgumentDiffFunction * > calculate_expected_unlocked_evolution(double phase_lag, double secondary_mass, bool decaying=true)
Calculate the predicted evolution for the test_unlocked_evolution() case.
Definition: testOrbitSolver.cpp:629
Core::AstroConst::solar_mass
const double solar_mass
Mass of the sun [kg].
Definition: AstronomicalConstants.h:21
Star::InterpolatedEvolutionStar::envelope
const EvolvingStellarEnvelope & envelope() const
The envelope of the star - inmodifiable.
Definition: EvolvingStar.h:94
AGE
AGE
The derivative w.r.t. age, excluding the dependence through the body&#39;s radius and the moments of iner...
Definition: DissipationQuantities.h:63
InCol::TSTART
The minimum age to start evolution at in Gyr.
Definition: IOColumns.h:67
PiecewiseFunction
Several functions stiched together.
Definition: PolynomialEvolution.h:477
Core::Error::General::get_message
const std::string & get_message()
Details about what caused the error.
Definition: Error.h:42
Evolve::BrokenPowerlawPhaseLagZone
A DissipatingZone where the phase lag is described by a broken powerlaw.
Definition: BrokenPowerlawPhaseLagZone.h:26
Evolve::ExpectedEvolutionMode::near_break
bool near_break(double age) const
Is the given age close to a break (hence ambigous mode).
Definition: ExpectedEvolutionMode.h:43
Core::EvolModeType
EvolModeType
The various evolution modes.
Definition: Common.h:42
Evolve::test_OrbitSolver::test_solution
void test_solution(const std::vector< const std::list< double > * > &tabulated_real_quantities, std::vector< const Core::OneArgumentDiffFunction *> expected_real_quantities, const ExpectedEvolutionMode< Core::EvolModeType > &expected_evol_mode, const ExpectedEvolutionMode< bool > &expected_wind_mode, double min_age, double max_age, bool debug_mode=false)
Tests the latest evolution calculated by the solver against the given tracks.
Definition: testOrbitSolver.cpp:303
Evolve::DissipatingZone::get_evolution_real
const std::list< double > & get_evolution_real(ZoneEvolutionQuantities quantity) const
The tabulated evolution of a real valued quantity so far.
Definition: DissipatingZone.h:835
Evolve::test_OrbitSolver::__temp_functions
std::vector< const Core::OneArgumentDiffFunction *> __temp_functions
A list of functions allocated during a test to delete at the end.
Definition: testOrbitSolver.h:51
Evolve::ANGULAR_MOMENTUM
Angular momentum of the zone.
Definition: DissipatingZone.h:30
Evolve::locked_unsat_eq_deriv
void locked_unsat_eq_deriv(double a, void *params, double *f, double *df)
The equation and its derivative that should be solved in order to get the semimamjor axis at a given ...
Definition: testOrbitSolver.cpp:1938
ExponentialPlusFunc
Represents a function of the form offset + scale*exp(rate*x) as well as its derivative.
Definition: PolynomialEvolution.h:394
Evolve::OrbitSolver
Solves the system of ODEs describing the evolution of a single planet around a single star...
Definition: OrbitSolver.h:115
Planet
Definition: Planet.h:15
Evolve::test_OrbitSolver::__solver
Evolve::OrbitSolver * __solver
The solver used for the current test.
Definition: testOrbitSolver.h:39
FunctionToPower
A function raised to some power.
Definition: PolynomialEvolution.h:556
Evolve::locked_unsat_eq
double locked_unsat_eq(double a, void *params)
The equation that should be solved in order to get the semimamjor axis at a gien time.
Definition: testOrbitSolver.cpp:1865
OutCol::CONV_PERIAPSIS
The orbital periapsis in the reference frame of the stellar convective zone in radians.
Definition: IOColumns.h:182
PiecewiseFunction::add_piece
void add_piece(const OneArgumentDiffFunction *piece)
Adds another piece of the function, where it applies is defined by its range.
Definition: PolynomialEvolution.cpp:471
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
FunctionRatio
The ratio of two functions;.
Definition: PolynomialEvolution.h:517
Evolve::test_OrbitSolver::get_evolution
std::vector< const std::list< double > * > get_evolution() const
Return the last calculated evolution.
Definition: testOrbitSolver.cpp:253
Evolve::locked_sat_deriv
double locked_sat_deriv(double a, void *params)
The derivative of the equation that should be solved in order to get the semimamjor axis at a gien ti...
Definition: testOrbitSolver.cpp:1925
Evolve::ExpectedEvolutionMode::add_break
void add_break(double age, MODE_TYPE mode)
Add an evolution mode that applies up to the given age.
Definition: ExpectedEvolutionMode.h:39
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::ExpectedEvolutionMode
Some evolution mode that changes at specified ages.
Definition: ExpectedEvolutionMode.h:22
Evolve::TransformedSolution
A class that can be passed to the solution testing function instead of the solver that transforms the...
Definition: TransformedSolution.h:32
Evolve::OrbitSolver::evolution_ages
const std::list< double > & evolution_ages() const
The ages at which evolution has been tabulated so far.
Definition: OrbitSolver.h:494
Evolve::TransformedSolution::add_transformation
void add_transformation(const std::vector< const Core::OneArgumentDiffFunction *> &transforms, double change_age)
Add more pieces to the transformation.
Definition: TransformedSolution.cpp:5
Core::AstroConst::G
const double G
Gravitational constant in SI.
Definition: AstronomicalConstants.h:26
Evolve::BrokenPowerlawPhaseLagZone::setup
void setup(const std::vector< double > &tidal_frequency_breaks, const std::vector< double > &spin_frequency_breaks, const std::vector< double > &tidal_frequency_powers, const std::vector< double > &spin_frequency_powers, double reference_phase_lag, double inertial_mode_enhancement=1.0, double inertial_mode_sharpness=10.0)
Seup the zone with the given breaks/powers and inertial mode enhancement. Continuous accress all brea...
Definition: BrokenPowerlawPhaseLagZone.cpp:218