Planetary Orbital Evolution due to Tides
Orbital evolution of two objects experiencing tides
DissipatingBody.cpp
1 #define BUILDING_LIBRARY
2 #include "DissipatingBody.h"
3 
4 namespace Evolve {
5 
6  double DissipatingBody::normalize_torques(double companion_mass,
7  double semimajor,
8  double orbital_frequency)
9  {
10  double torque_norm = (
11  std::pow(companion_mass / std::pow(semimajor, 3), 2)
12  *
13  std::pow(radius(), 5)
14  *
15  Core::AstroConst::G
16  *
17  Core::AstroConst::day
18  *
19  Core::AstroConst::Gyr
20  *
21  Core::AstroConst::solar_mass
22  /
23  std::pow(Core::AstroConst::solar_radius, 3)
24  );
25  double dangular_velocity_da = -1.5 * orbital_frequency / semimajor;
26  for(unsigned zone_index = 0; zone_index < number_zones(); ++zone_index) {
27  bool above = false;
28  const DissipatingZone &current_zone = zone(zone_index);
29  do {
30  std::valarray<Eigen::Vector3d> &tidal_torque = (
31  above
32  ? __tidal_torques_above
33  : __tidal_torques_below
34  )[zone_index];
35  tidal_torque[Dissipation::ORBITAL_FREQUENCY] += (
36  4.0
37  /
38  orbital_frequency
39  *
40  tidal_torque[Dissipation::NO_DERIV]
41  );
42  tidal_torque[Dissipation::RADIUS] = (
43  5.0 / radius() * tidal_torque[Dissipation::NO_DERIV]
44  );
45  tidal_torque[Dissipation::MOMENT_OF_INERTIA] = (
46  -current_zone.spin_frequency()
47  /
48  current_zone.moment_of_inertia()
49  *
50  tidal_torque[Dissipation::SPIN_FREQUENCY]
51  );
52  tidal_torque[Dissipation::SPIN_ANGMOM] = (
53  tidal_torque[Dissipation::SPIN_FREQUENCY]
54  /
55  current_zone.moment_of_inertia()
56  );
57  tidal_torque[Dissipation::SEMIMAJOR] = (
58  dangular_velocity_da
59  *
60  tidal_torque[Dissipation::ORBITAL_FREQUENCY]
61  );
62  above = !above;
63  } while (above);
64  }
65 
66  assert(__tidal_torques_above.size() == __tidal_torques_below.size());
67  for(unsigned i = 0; i < __tidal_torques_above.size(); ++i) {
68  assert(__tidal_torques_above[i].size()
69  ==
70  __tidal_torques_below[i].size());
71  for(unsigned j = 0; j < __tidal_torques_above[i].size(); ++j) {
72  __tidal_torques_above[i][j] *= torque_norm;
73  __tidal_torques_below[i][j] *= torque_norm;
74  }
75  }
76  return torque_norm;
77  }
78 
79  void DissipatingBody::collect_orbit_rates(double orbital_frequency,
80  double torque_norm)
81  {
82  __power_norm = torque_norm * orbital_frequency;
83  DissipatingZone &surface_zone = zone(0);
84  unsigned invalid_ind = __orbit_entries.size(),
85  no_deriv_ind = invalid_ind,
86  orbital_freq_deriv_ind = invalid_ind,
87  radius_deriv_ind = invalid_ind,
88  semimajor_deriv_ind = invalid_ind;
89  for(unsigned i=0; i<__orbit_entries.size(); ++i)
90  switch(__orbit_entries[i]) {
91  case Dissipation::NO_DERIV : no_deriv_ind=i; break;
92  case Dissipation::ORBITAL_FREQUENCY : orbital_freq_deriv_ind=i;
93  break;
94  case Dissipation::RADIUS : radius_deriv_ind=i; break;
95  case Dissipation::SEMIMAJOR : semimajor_deriv_ind=i; break;
96  default:;
97  }
98  __orbit_power = 0;
99  for(
100  unsigned zone_index = 0;
101  zone_index < number_zones();
102  ++zone_index
103  ) {
104  DissipatingZone &current_zone = zone(zone_index);
105  std::valarray<Eigen::Vector3d> &tidal_torque = (
106  __tidal_torques_below[zone_index]
107  );
108  for(
109  unsigned deriv_ind = 0;
110  deriv_ind < __orbit_entries.size();
111  ++deriv_ind
112  ) {
113  Dissipation::QuantityEntry entry = __orbit_entries[deriv_ind];
114  if(entry == Dissipation::EXPANSION_ERROR) {
115  __orbit_power[deriv_ind] += std::abs(
116  current_zone.tidal_power(false,
117  Dissipation::EXPANSION_ERROR)
118  );
119  __orbit_torque[entry].setConstant(tidal_torque[entry].norm());
120  }
121  else if(entry < Dissipation::END_DIMENSIONLESS_DERIV)
122  __orbit_power[deriv_ind] -= (
123  current_zone.tidal_power(false, entry)
124  );
125  if(zone_index) {
126  __orbit_torque[entry] -= zone_to_zone_transform(
127  current_zone,
128  surface_zone,
129  tidal_torque[entry]
130  );
131  } else __orbit_torque[entry] = -tidal_torque[entry];
132  }
133  if(zone_index) {
134  __orbit_torque[Dissipation::INCLINATION] -= (
135  zone_to_zone_transform(current_zone,
136  surface_zone,
137  tidal_torque[Dissipation::NO_DERIV],
138  Dissipation::INCLINATION,
139  false)
140  );
141  __orbit_torque[Dissipation::PERIAPSIS] -= (
142  zone_to_zone_transform(current_zone,
143  surface_zone,
144  tidal_torque[Dissipation::NO_DERIV],
145  Dissipation::PERIAPSIS,
146  false)
147  );
148  }
149  }
150  __orbit_power *= __power_norm;
151  __orbit_power[orbital_freq_deriv_ind] += (
152  5.0 / orbital_frequency * __orbit_power[no_deriv_ind]
153  );
154  __orbit_power[radius_deriv_ind] += (
155  5.0 / radius() * __orbit_power[no_deriv_ind]
156  );
157  __orbit_power[semimajor_deriv_ind] = (
158  __dorbital_frequency_da
159  *
160  __orbit_power[orbital_freq_deriv_ind]
161  );
162  }
163 
164  void DissipatingBody::calculate_orbit_rate_corrections()
165  {
166  __orbit_power_correction.resize(__num_locked_zones);
167  __orbit_torque_correction.resize(__num_locked_zones);
168  unsigned correction_index = 0;
169  for(unsigned zone_index=0; zone_index<number_zones(); ++zone_index)
170  if(zone(zone_index).locked()) {
171  DissipatingZone &this_zone=zone(zone_index);
172  __orbit_power_correction[correction_index] = (
173  tidal_power(zone_index, false)
174  -
175  tidal_power(zone_index, true)
176  );
177  __orbit_torque_correction[correction_index] = (
178  zone_to_zone_transform(this_zone, zone(0),
179  tidal_torque(zone_index, false)
180  -
181  tidal_torque(zone_index, true))
182  );
183  ++correction_index;
184  }
185  }
186 
187  void DissipatingBody::angular_momentum_transfer(
188  const DissipatingZone &outer_zone,
189  const DissipatingZone &inner_zone,
190  Eigen::Vector3d &outer_angmom_gain,
191  Eigen::Vector3d &inner_angmom_gain,
192  Dissipation::QuantityEntry deriv,
193  bool with_respect_to_outer
194  ) const
195  {
196  assert(deriv == Dissipation::NO_DERIV
197  ||
198  deriv == Dissipation::INCLINATION
199  ||
200  deriv == Dissipation::PERIAPSIS);
201 
202  double dm_dt = inner_zone.outer_mass(1),
203  lost_spin = (dm_dt >= 0
204  ? outer_zone
205  : inner_zone).spin_frequency(),
206  angmom_transfer = (
207  -2.0 / 3.0
208  *
209  std::pow(inner_zone.outer_radius(), 2)
210  *
211  lost_spin * dm_dt
212  );
213  if(dm_dt > 0) {
214  outer_angmom_gain[0] = outer_angmom_gain[1] = 0;
215  outer_angmom_gain[2] = (deriv == Dissipation::NO_DERIV
216  ? angmom_transfer
217  : 0);
218  inner_angmom_gain = -zone_to_zone_transform(
219  outer_zone,
220  inner_zone,
221  Eigen::Vector3d(0, 0, angmom_transfer),
222  deriv,
223  with_respect_to_outer
224  );
225  } else {
226  inner_angmom_gain[0] = inner_angmom_gain[1] = 0;
227  inner_angmom_gain[2] = (deriv == Dissipation::NO_DERIV
228  ? -angmom_transfer
229  : 0);
230  outer_angmom_gain=zone_to_zone_transform(
231  inner_zone,
232  outer_zone,
233  Eigen::Vector3d(0, 0, angmom_transfer),
234  deriv,
235  !with_respect_to_outer
236  );
237  }
238  }
239 
240  Eigen::Vector3d DissipatingBody::angular_momentum_transfer_from_top(
241  unsigned zone_index,
242  Dissipation::QuantityEntry deriv,
243  bool with_respect_to_outer
244  ) const
245  {
246  assert(zone_index > 0);
247 
248  const DissipatingZone &this_zone = zone(zone_index),
249  &zone_above = zone(zone_index - 1);
250  double scaling = Core::NaN,
251  dm_dt = this_zone.outer_mass(1);
252  if(deriv == Dissipation::NO_DERIV)
253  scaling = 1;
254  else if(deriv == Dissipation::AGE) {
255  scaling = (
256  2.0 * this_zone.outer_radius(1) / this_zone.outer_radius(0)
257  +
258  this_zone.outer_mass(2) / dm_dt
259  );
260  } else if(
261  deriv == Dissipation::SPIN_FREQUENCY
262  ||
263  deriv == Dissipation::MOMENT_OF_INERTIA
264  ||
265  deriv == Dissipation::SPIN_ANGMOM
266  ) {
267  if(
268  (with_respect_to_outer && dm_dt <= 0)
269  ||
270  (!with_respect_to_outer && dm_dt>=0)
271  )
272  return Eigen::Vector3d(0, 0, 0);
273  else if(deriv == Dissipation::SPIN_FREQUENCY)
274  scaling = 1.0 / (with_respect_to_outer
275  ? zone_above.spin_frequency()
276  : this_zone.spin_frequency());
277  else if(deriv == Dissipation::MOMENT_OF_INERTIA)
278  scaling = -1.0 / (with_respect_to_outer
279  ? zone_above.moment_of_inertia()
280  : this_zone.moment_of_inertia());
281  else scaling = 1.0 / (with_respect_to_outer
282  ? zone_above.angular_momentum()
283  : this_zone.angular_momentum());
284  } else if(
285  deriv == Dissipation::INCLINATION
286  ||
287  deriv == Dissipation::PERIAPSIS
288  ) {
289  if(dm_dt <= 0)
290  return Eigen::Vector3d(0, 0, 0);
291  else {
292  Eigen::Vector3d dummy, result;
293  angular_momentum_transfer(zone_above,
294  this_zone,
295  dummy,
296  result,
297  deriv,
298  with_respect_to_outer);
299  return result;
300  }
301  } else
302  assert(false);
303  return scaling * __angular_momentum_transfer[zone_index - 1][1];
304  }
305 
306  Eigen::Vector3d DissipatingBody::angular_momentum_transfer_from_bottom(
307  unsigned zone_index,
308  Dissipation::QuantityEntry deriv,
309  bool with_respect_to_inner
310  ) const
311  {
312  assert(zone_index < number_zones() - 1);
313 
314  const DissipatingZone &this_zone = zone(zone_index),
315  &zone_below = zone(zone_index + 1);
316  double scaling = Core::NaN,
317  dm_dt = zone_below.outer_mass(1);
318  if(deriv == Dissipation::NO_DERIV)
319  scaling = 1;
320  else if(deriv == Dissipation::AGE) {
321  scaling = (
322  2.0 * zone_below.outer_radius(1) / zone_below.outer_radius(0)
323  +
324  zone_below.outer_mass(2) / dm_dt
325  );
326  } else if(
327  deriv == Dissipation::SPIN_FREQUENCY
328  ||
329  deriv == Dissipation::MOMENT_OF_INERTIA
330  ||
331  deriv == Dissipation::SPIN_ANGMOM
332  ) {
333  if(
334  (with_respect_to_inner && dm_dt >= 0)
335  ||
336  (!with_respect_to_inner && dm_dt <= 0)
337  )
338  return Eigen::Vector3d(0, 0, 0);
339  else if(deriv == Dissipation::SPIN_FREQUENCY)
340  scaling = 1.0 / (with_respect_to_inner
341  ? zone_below.spin_frequency()
342  : this_zone.spin_frequency());
343  else if(deriv == Dissipation::MOMENT_OF_INERTIA)
344  scaling = -1.0 / (with_respect_to_inner
345  ? zone_below.moment_of_inertia()
346  : this_zone.moment_of_inertia());
347  else scaling = 1.0 / (with_respect_to_inner
348  ? zone_below.angular_momentum()
349  : this_zone.angular_momentum());
350  } else if(
351  deriv == Dissipation::INCLINATION
352  ||
353  deriv == Dissipation::PERIAPSIS
354  ) {
355  if(dm_dt >= 0)
356  return Eigen::Vector3d(0, 0, 0);
357  else {
358  Eigen::Vector3d dummy, result;
359  angular_momentum_transfer(this_zone,
360  zone_below,
361  result,
362  dummy,
363  deriv,
364  !with_respect_to_inner);
365  }
366  } else
367  assert(false);
368 
369  return scaling * __angular_momentum_transfer[zone_index][0];
370  }
371 
372  Eigen::Vector3d DissipatingBody::angular_momentum_transfer_to_zone(
373  unsigned zone_index,
374  Dissipation::QuantityEntry deriv,
375  int deriv_zone
376  ) const
377  {
378 
379  assert(deriv < Dissipation::NUM_DERIVATIVES);
380  if(
381  deriv == Dissipation::ORBITAL_FREQUENCY
382  ||
383  deriv == Dissipation::ECCENTRICITY
384  ||
385  deriv == Dissipation::RADIUS
386  ||
387  deriv == Dissipation::SEMIMAJOR
388  )
389  return Eigen::Vector3d(0,0,0);
390  Eigen::Vector3d result(0, 0, 0);
391  if(
392  zone_index > 0
393  &&
394  (
395  deriv == Dissipation::NO_DERIV
396  ||
397  deriv == Dissipation::AGE
398  ||
399  deriv_zone <= 0
400  )
401  )
402  result = angular_momentum_transfer_from_top(zone_index,
403  deriv,
404  deriv_zone < 0);
405  if(
406  zone_index < number_zones() - 1
407  &&
408  (
409  deriv == Dissipation::NO_DERIV
410  ||
411  deriv == Dissipation::AGE
412  ||
413  deriv_zone >= 0
414  )
415  )
416  result += angular_momentum_transfer_from_bottom(zone_index,
417  deriv,
418  deriv_zone > 0);
419  return result;
420  }
421 
422  DissipatingBody::DissipatingBody() :
423  __orbit_entries(7),
424  __orbit_power(__orbit_entries.size()),
425  __orbit_torque(Dissipation::NUM_ENTRIES),
426  __orbit_torque_correction(0),
427  __num_locked_zones(0)
428  {
429  __orbit_entries[0] = Dissipation::NO_DERIV;
430  __orbit_entries[1] = Dissipation::AGE;
431  __orbit_entries[2] = Dissipation::ORBITAL_FREQUENCY;
432  __orbit_entries[3] = Dissipation::ECCENTRICITY;
433  __orbit_entries[4] = Dissipation::RADIUS;
434  __orbit_entries[5] = Dissipation::SEMIMAJOR;
435  __orbit_entries[6] = Dissipation::EXPANSION_ERROR;
436  }
437 
438  void DissipatingBody::correct_orbit_power(
439  Eigen::VectorXd &above_lock_fractions_age_deriv,
440  Eigen::VectorXd &above_lock_fractions_semimajor_deriv,
441  Eigen::VectorXd &above_lock_fractions_eccentricity_deriv,
442  Eigen::VectorXd &above_lock_fractions_radius_deriv
443  )
444  {
445  unsigned correction_index = 0;
446  for(
447  unsigned zone_index = 0;
448  zone_index < number_zones();
449  ++zone_index
450  ) {
451  if(zone(zone_index).locked()) {
452  unsigned locked_zone_index = (
453  zone(zone_index).locked_zone_index()
454  );
455  __orbit_power_correction[correction_index] = (
456  tidal_power(zone_index, false)
457  -
458  tidal_power(zone_index, true)
459  );
460  for(
461  unsigned entry_ind = 0;
462  entry_ind < __orbit_entries.size();
463  ++entry_ind
464  ) {
465  Dissipation::QuantityEntry entry = (
466  __orbit_entries[entry_ind]
467  );
468  __orbit_power[entry_ind] += (
469  __above_lock_fractions
470  [Dissipation::NO_DERIV]
471  [locked_zone_index]
472  *
473  (
474  tidal_power(zone_index, false, entry)
475  -
476  tidal_power(zone_index, true, entry)
477  )
478  );
479  if(
480  entry != Dissipation::NO_DERIV
481  &&
482  entry != Dissipation::EXPANSION_ERROR
483  ) {
484  double frac_deriv = Core::NaN;
485  if(entry == Dissipation::AGE)
486  frac_deriv = (above_lock_fractions_age_deriv
487  [locked_zone_index]);
488  else if(entry == Dissipation::ORBITAL_FREQUENCY)
489  frac_deriv = (
490  above_lock_fractions_semimajor_deriv
491  [locked_zone_index]
492  /
493  __dorbital_frequency_da
494  );
495  else if(entry == Dissipation::ECCENTRICITY)
496  frac_deriv = (
497  above_lock_fractions_eccentricity_deriv
498  [locked_zone_index]
499  /
500  __dorbital_frequency_da
501  );
502  else if(entry == Dissipation::RADIUS)
503  frac_deriv = (above_lock_fractions_radius_deriv
504  [locked_zone_index]);
505  else if(entry == Dissipation::SEMIMAJOR)
506  frac_deriv = (above_lock_fractions_semimajor_deriv
507  [locked_zone_index]);
508  else
509  assert(false);
510 
511  __orbit_power[entry_ind] += (
512  frac_deriv
513  *
514  __orbit_power_correction[correction_index]
515  );
516  }
517  }
518  ++correction_index;
519  }
520  }
521  }
522 
523  void DissipatingBody::correct_orbit_torque(
524  std::valarray<Eigen::VectorXd> &above_lock_fractions
525  )
526  {
527  DissipatingZone &surface_zone = zone(0);
528  unsigned correction_index = 0;
529  for(
530  unsigned zone_index = 0;
531  zone_index < number_zones();
532  ++zone_index
533  ) {
534  if(zone(zone_index).locked()) {
535  unsigned locked_zone_index = (
536  zone(zone_index).locked_zone_index()
537  );
538  DissipatingZone &this_zone = zone(zone_index);
539  for(
540  int entry_int = Dissipation::NO_DERIV;
541  entry_int < Dissipation::NUM_ENTRIES;
542  ++entry_int
543  ) {
544  Dissipation::QuantityEntry entry = (
545  static_cast<Dissipation::QuantityEntry>(entry_int)
546  );
547  Eigen::Vector3d correction = (
548  above_lock_fractions
549  [Dissipation::NO_DERIV]
550  [locked_zone_index]
551  *
552  (
553  tidal_torque(zone_index, false, entry)
554  -
555  tidal_torque(zone_index, true, entry)
556  )
557  );
558  if(entry == Dissipation::EXPANSION_ERROR) {
559  __orbit_torque[entry].array() += correction.norm();
560  } else {
561  if(zone_index) {
562  __orbit_torque[entry] += zone_to_zone_transform(
563  this_zone,
564  surface_zone,
565  correction
566  );
567  if(
568  entry == Dissipation::INCLINATION
569  ||
570  entry == Dissipation::PERIAPSIS
571  )
572  __orbit_torque[entry] += zone_to_zone_transform(
573  this_zone,
574  surface_zone,
575  (
576  above_lock_fractions
577  [Dissipation::NO_DERIV]
578  [locked_zone_index]
579  *
580  __orbit_torque_correction
581  [correction_index]
582  ),
583  static_cast<Dissipation::QuantityEntry>(
584  entry
585  )
586  );
587  } else
588  __orbit_torque[entry] += correction;
589 
590  if(entry != Dissipation::NO_DERIV)
591  __orbit_torque[entry] += (
592  above_lock_fractions[entry][locked_zone_index]
593  *
594  __orbit_torque_correction[correction_index]
595  );
596  }
597  }
598  ++correction_index;
599  }
600  }
601  }
602 
603  void DissipatingBody::configure(bool initialize,
604  double age,
605  double companion_mass,
606  double semimajor,
607  double eccentricity,
608  const double *spin_angmom,
609  const double *inclination,
610  const double *periapsis,
611  bool locked_surface,
612  bool zero_outer_inclination,
613  bool zero_outer_periapsis)
614  {
615  double orbital_angmom=Core::orbital_angular_momentum(
616  companion_mass,
617  mass(),
618  semimajor,
619  eccentricity
620  );
621  __orbital_frequency=Core::orbital_angular_velocity(
622  companion_mass,
623  mass(),
624  semimajor,
625  false
626  );
627  __dorbital_frequency_da=Core::orbital_angular_velocity(
628  companion_mass,
629  mass(),
630  semimajor,
631  true
632  );
633 
634  if(initialize) {
635  __num_locked_zones = 0;
636 #ifndef NDEBUG
637  std::cerr << "Initializing DissipatingBody" << std::endl;
638 #endif
639  }
640 
641  __tidal_torques_above.resize(number_zones());
642  __tidal_torques_below.resize(number_zones());
643  __angular_momentum_transfer.resize(number_zones() - 1);
644  unsigned angmom_offset = (locked_surface ? 1 : 0);
645  for(
646  unsigned zone_index = 0;
647  zone_index < number_zones();
648  ++zone_index
649  ) {
650  DissipatingZone &current_zone = zone(zone_index);
651  double zone_inclination, zone_periapsis, zone_spin;
652  if(!inclination)
653  zone_inclination = 0;
654  else if(zero_outer_inclination)
655  zone_inclination = (zone_index
656  ? inclination[zone_index - 1]
657  : 0);
658  else
659  zone_inclination = inclination[zone_index];
660  if(!periapsis)
661  zone_periapsis = 0;
662  else if(zero_outer_periapsis)
663  zone_periapsis = (zone_index ? periapsis[zone_index-1] : 0);
664  else
665  zone_periapsis = periapsis[zone_index];
666  if(locked_surface && zone_index == 0)
667  zone_spin = surface_lock_frequency();
668  else if(current_zone.locked() && !initialize) {
669  zone_spin = Core::NaN;
670  ++angmom_offset;
671  } else
672  zone_spin = spin_angmom[zone_index - angmom_offset];
673 #ifndef NDEBUG
674  std::cerr << "Configuring zone " << zone_index << std::endl;
675 #endif
676  current_zone.configure(initialize,
677  age,
678  __orbital_frequency,
679  eccentricity,
680  orbital_angmom,
681  zone_spin,
682  zone_inclination,
683  zone_periapsis,
684  locked_surface && zone_index == 0);
685  }
686  for(
687  unsigned zone_index = 0;
688  zone_index < number_zones();
689  ++zone_index
690  ) {
691  DissipatingZone &current_zone = zone(zone_index);
692  if(zone_index < number_zones() - 1) {
693  __angular_momentum_transfer[zone_index].resize(2);
694  angular_momentum_transfer(
695  current_zone,
696  zone(zone_index + 1),
697  __angular_momentum_transfer[zone_index][0],
698  __angular_momentum_transfer[zone_index][1]
699  );
700  }
701  bool above = false;
702  do {
703  std::valarray<Eigen::Vector3d> &tidal_torque =
704  (above ? __tidal_torques_above : __tidal_torques_below)
705  [zone_index];
706  tidal_torque.resize(Dissipation::NUM_ENTRIES);
707  for(
708  int torque_ind = Dissipation::NO_DERIV;
709  torque_ind <= Dissipation::END_DIMENSIONLESS_DERIV;
710  ++torque_ind
711  ) {
712  Dissipation::QuantityEntry entry = (
713  torque_ind < Dissipation::END_DIMENSIONLESS_DERIV
714  ? static_cast<Dissipation::QuantityEntry>(torque_ind)
715  : Dissipation::EXPANSION_ERROR
716  );
717  tidal_torque[entry][0] = current_zone.tidal_torque_x(
718  above,
719  entry
720  );
721  tidal_torque[entry][1] = current_zone.tidal_torque_y(
722  above,
723  entry
724  );
725  tidal_torque[entry][2] = current_zone.tidal_torque_z(
726  above,
727  entry
728  );
729  assert(!std::isnan(tidal_torque[entry].sum()));
730  }
731  above = !above;
732  } while(above);
733  }
734  collect_orbit_rates(__orbital_frequency,
735  normalize_torques(companion_mass,
736  semimajor,
737  __orbital_frequency));
738  calculate_orbit_rate_corrections();
739  __above_lock_fractions.resize(0);
740  }
741 
742  Eigen::Vector3d DissipatingBody::nontidal_torque(
743  unsigned zone_index,
744  Dissipation::QuantityEntry deriv,
745  int deriv_zone
746  ) const
747  {
748  assert(deriv < Dissipation::NUM_DERIVATIVES);
749  assert(zone_index < number_zones());
750  assert(static_cast<int>(zone_index) + deriv_zone >= 0);
751  assert(static_cast<int>(zone_index) + deriv_zone
752  <
753  static_cast<int>(number_zones()));
754  const DissipatingZone &this_zone = zone(zone_index);
755  Eigen::Vector3d result(0, 0, 0);
756  if(zone_index == 0 && deriv_zone == 0)
757  result[2] = -angular_momentum_loss(deriv);
758  result += angular_momentum_transfer_to_zone(zone_index,
759  deriv,
760  deriv_zone);
761  if(
762  zone_index < number_zones() - 1
763  &&
764  (!zone_specific(deriv) || deriv_zone >= 0)
765  )
766  result += angular_momentum_coupling(zone_index,
767  deriv,
768  deriv_zone == 0);
769  if(
770  zone_index > 0
771  &&
772  (!zone_specific(deriv) || deriv_zone<=0)
773  ) {
774  result -= zone_to_zone_transform(
775  zone(zone_index-1),
776  this_zone,
777  angular_momentum_coupling(zone_index - 1,
778  deriv,
779  deriv_zone < 0)
780  );
781  if(
782  deriv == Dissipation::INCLINATION
783  ||
784  deriv == Dissipation::PERIAPSIS
785  )
786  result -= zone_to_zone_transform(
787  zone(zone_index - 1),
788  this_zone,
789  angular_momentum_coupling(zone_index - 1),
790  deriv,
791  deriv_zone<0
792  );
793  }
794  return result;
795  }
796 
797  double DissipatingBody::tidal_power(unsigned zone_index,
798  bool above,
799  Dissipation::QuantityEntry entry) const
800  {
801  assert(zone_index < number_zones());
802 
803  const DissipatingZone &this_zone = zone(zone_index);
804  double result = (entry < Dissipation::END_DIMENSIONLESS_DERIV
805  ? __power_norm * this_zone.tidal_power(above, entry)
806  : 0.0);
807  if(
808  entry == Dissipation::ORBITAL_FREQUENCY
809  ||
810  entry == Dissipation::SEMIMAJOR
811  ) {
812  result += (5.0
813  /
814  __orbital_frequency
815  *
816  __power_norm
817  *
818  this_zone.tidal_power(above));
819  if(entry == Dissipation::SEMIMAJOR)
820  result *= __dorbital_frequency_da;
821  }
822  else if(entry == Dissipation::RADIUS)
823  result += 5.0 / radius() * this_zone.tidal_power(above);
824  return result;
825  }
826 
827  void DissipatingBody::set_above_lock_fractions(
828  std::valarray<Eigen::VectorXd> &above_lock_fractions
829  )
830  {
831  assert(above_lock_fractions.size() == Dissipation::NUM_DERIVATIVES);
832 
833  __above_lock_fractions.resize(Dissipation::NUM_DERIVATIVES);
834  for(unsigned i = 0; i < Dissipation::NUM_DERIVATIVES; ++i) {
835  assert(above_lock_fractions[i].size() >= __num_locked_zones);
836 
837  __above_lock_fractions[i].resize(above_lock_fractions[i].size());
838  }
839  __above_lock_fractions = above_lock_fractions;
840  correct_orbit_power(
841  above_lock_fractions[Dissipation::AGE],
842  above_lock_fractions[Dissipation::SEMIMAJOR],
843  above_lock_fractions[Dissipation::ECCENTRICITY],
844  above_lock_fractions[Dissipation::RADIUS]
845  );
846  correct_orbit_torque(above_lock_fractions);
847  }
848 
849  double DissipatingBody::tidal_orbit_power(
850  Dissipation::QuantityEntry entry,
851  unsigned deriv_zone_index,
852  const Eigen::VectorXd &above_lock_fraction_deriv
853  ) const
854  {
855  const DissipatingZone &deriv_zone = zone(deriv_zone_index);
856  double result = Core::NaN;
857  if(
858  !zone_specific(entry)
859  ||
860  !deriv_zone.locked()
861  ||
862  __above_lock_fractions.size() == 0
863  ) {
864  if(zone_specific(entry)) {
865  result = tidal_power(deriv_zone_index, false, entry);
866  } else {
867  for(unsigned i = 0; i < __orbit_entries.size(); ++i)
868  if(entry == __orbit_entries[i])
869  return __orbit_power[i];
870 
871  assert(false);
872  }
873  if(__above_lock_fractions.size() == 0)
874  return result;
875  } else {
876  double above_frac = (__above_lock_fractions
877  [Dissipation::NO_DERIV]
878  [zone(deriv_zone_index).locked_zone_index()]);
879  result = (above_frac
880  *
881  tidal_power(deriv_zone_index, true, entry)
882  +
883  (1.0 - above_frac)
884  *
885  tidal_power(deriv_zone_index, false, entry));
886  }
887 
888  assert(above_lock_fraction_deriv.size()
889  ==
890  __above_lock_fractions[Dissipation::NO_DERIV].size());
891 
892  unsigned correction_index = 0;
893  for(unsigned zone_index = 0; zone_index < number_zones(); ++zone_index)
894  if(zone(zone_index).locked()) {
895  result += (
896  above_lock_fraction_deriv
897  [zone(zone_index).locked_zone_index()]
898  *
899  __orbit_power_correction[correction_index]
900  );
901  ++correction_index;
902  }
903  return result;
904  }
905 
906  Eigen::Vector3d DissipatingBody::tidal_orbit_torque(
907  Dissipation::QuantityEntry entry,
908  unsigned deriv_zone_index,
909  const Eigen::VectorXd &above_lock_fraction_deriv
910  ) const
911  {
912  if(zone_specific(entry) && deriv_zone_index != 0) {
913  const DissipatingZone &deriv_zone = zone(deriv_zone_index);
914  double above_frac = (__above_lock_fractions
915  [Dissipation::NO_DERIV]
916  [deriv_zone.locked_zone_index()]);
917  Eigen::Vector3d result = zone_to_zone_transform(
918  deriv_zone,
919  zone(0),
920  (
921  (above_frac - 1.0)
922  *
923  __tidal_torques_below[deriv_zone_index][entry]
924  -
925  above_frac
926  *
927  __tidal_torques_above[deriv_zone_index][entry]
928  )
929  );
930 
931  assert(above_lock_fraction_deriv.size()
932  ==
933  __above_lock_fractions[Dissipation::NO_DERIV].size());
934 
935  unsigned correction_index = 0;
936  for(
937  unsigned zone_index = 0;
938  zone_index < number_zones();
939  ++zone_index
940  )
941  if(zone(zone_index).locked()) {
942  result += (
943  above_lock_fraction_deriv[
944  zone(zone_index).locked_zone_index()
945  ]
946  *
947  __orbit_torque_correction[correction_index]
948  );
949  ++correction_index;
950  }
951  if(
952  entry == Dissipation::INCLINATION
953  ||
954  entry == Dissipation::PERIAPSIS
955  )
956  result += zone_to_zone_transform(
957  deriv_zone,
958  zone(0),
959  (
960  above_frac
961  *
962  (
963  __tidal_torques_above
964  [deriv_zone_index]
965  [Dissipation::NO_DERIV]
966  )
967  +
968  (1.0 - above_frac)
969  *
970  (
971  __tidal_torques_below
972  [deriv_zone_index]
973  [Dissipation::NO_DERIV]
974  )
975  ),
976  entry,
977  true
978  );
979  return result;
980  } else {
981  return __orbit_torque[entry];
982  }
983  }
984 
985  Eigen::Vector3d DissipatingBody::tidal_orbit_torque(
986  const DissipatingZone &reference_zone,
987  Dissipation::QuantityEntry entry,
988  unsigned deriv_zone_index,
989  const Eigen::VectorXd &above_lock_fraction_deriv
990  ) const
991  {
992  Eigen::Vector3d result;
993  if(entry == Dissipation::EXPANSION_ERROR) {
994  result = tidal_orbit_torque(Dissipation::EXPANSION_ERROR,
995  0,
996  Eigen::VectorXd());
997  } else {
998  result = zone_to_zone_transform(
999  zone(0),
1000  reference_zone,
1001  tidal_orbit_torque(entry,
1002  deriv_zone_index,
1003  above_lock_fraction_deriv)
1004  );
1005  if(
1006  (
1007  entry == Dissipation::INCLINATION
1008  ||
1009  entry == Dissipation::PERIAPSIS
1010  )
1011  &&
1012  (
1013  deriv_zone_index == 0
1014  ||
1015  &zone(deriv_zone_index) == &reference_zone
1016  )
1017  ) {
1018  result += zone_to_zone_transform(zone(0),
1019  reference_zone,
1020  tidal_orbit_torque(),
1021  entry,
1022  deriv_zone_index==0);
1023  }
1024  }
1025  return result;
1026  }
1027 
1028  void DissipatingBody::add_to_evolution()
1029  {
1030  for(unsigned zone_ind = 0; zone_ind < number_zones(); ++zone_ind)
1031  zone(zone_ind).add_to_evolution();
1032  }
1033 
1034  void DissipatingBody::rewind_evolution(unsigned nsteps)
1035  {
1036  for(unsigned zone_ind = 0; zone_ind < number_zones(); ++zone_ind)
1037  zone(zone_ind).rewind_evolution(nsteps);
1038  }
1039 
1040  void DissipatingBody::reset_evolution()
1041  {
1042  for(unsigned zone_ind = 0; zone_ind < number_zones(); ++zone_ind)
1043  zone(zone_ind).reset_evolution();
1044  }
1045 
1046  CombinedStoppingCondition *DissipatingBody::stopping_conditions(
1047  BinarySystem &system,
1048  bool primary
1049  )
1050  {
1051  CombinedStoppingCondition *result = new CombinedStoppingCondition();
1052  for(unsigned zone_ind = 0; zone_ind < number_zones(); ++zone_ind)
1053  (*result) |= zone(zone_ind).stopping_conditions(system,
1054  primary,
1055  zone_ind);
1056  return result;
1057  }
1058 
1059  void DissipatingBody::change_e_order(unsigned new_e_order,
1060  BinarySystem &system,
1061  bool primary)
1062  {
1063  for(unsigned zone_ind = 0; zone_ind < number_zones(); ++zone_ind)
1064  if(zone(zone_ind).dissipative())
1065  zone(zone_ind).change_e_order(new_e_order,
1066  system,
1067  primary,
1068  zone_ind);
1069  }
1070 
1071 } //End Envolve namespace.
Evolve::DissipatingBody::__power_norm
double __power_norm
The coefficient used to normalize tidal power.
Definition: DissipatingBody.h:43
Evolve::DissipatingZone::tidal_power
double tidal_power(bool above, Dissipation::QuantityEntry entry=Dissipation::NO_DERIV) const
The dimensionless tidal power or one of its derivatives.
Definition: DissipatingZone.h:590
Evolve::DissipatingZone::locked_zone_index
unsigned locked_zone_index() const
The index of this zone in the list of locked zones (valid only if locked).
Definition: DissipatingZone.h:857
Evolve::DissipatingBody::__angular_momentum_transfer
std::valarray< std::valarray< Eigen::Vector3d > > __angular_momentum_transfer
The rate of angular momentum transfer between two neighboring zones due to zone boundaries moving...
Definition: DissipatingBody.h:63
Evolve::DissipatingBody::rewind_evolution
virtual void rewind_evolution(unsigned nsteps)
Discards the last steps from the evolution.
Definition: DissipatingBody.cpp:1034
Evolve::DissipatingBody::reset_evolution
virtual void reset_evolution()
Discards all evolution.
Definition: DissipatingBody.cpp:1040
RADIUS
RADIUS
The derivative w.r.t. the radius of the body in .
Definition: DissipationQuantities.h:99
Evolve::DissipatingZone::configure
virtual void configure(bool initialize, double age, double orbital_frequency, double eccentricity, double orbital_angmom, double spin, double inclination, double periapsis, bool spin_is_frequency)
Defines the current orbit, triggering re-calculation of all quantities.
Definition: DissipatingZone.cpp:485
Evolve::DissipatingBody::add_to_evolution
virtual void add_to_evolution()
Appends the state defined by last configure(), to the evolution.
Definition: DissipatingBody.cpp:1028
Evolve::DissipatingZone::angular_momentum
double angular_momentum() const
The angular momentum of the given zone in .
Definition: DissipatingZone.h:587
Evolve::DissipatingBody::set_above_lock_fractions
void set_above_lock_fractions(std::valarray< Eigen::VectorXd > &above_lock_fractions)
Corrects the tidal orbit energy gain and angular momentum gain for locked zones.
Definition: DissipatingBody.cpp:827
Evolve::DissipatingBody::angular_momentum_transfer_from_bottom
Eigen::Vector3d angular_momentum_transfer_from_bottom(unsigned zone_index, Dissipation::QuantityEntry deriv=Dissipation::NO_DERIV, bool with_respect_to_inner=false) const
Rate of angular momentum transfer (or its derivatives) to a zone due to its bottom boundary moving...
Definition: DissipatingBody.cpp:306
SPIN_FREQUENCY
SPIN_FREQUENCY
The derivative w.r.t. the spin frequency of a dissipating zone.
Definition: DissipationQuantities.h:63
Evolve::DissipatingBody::DissipatingBody
DissipatingBody()
Some initializations for new objects.
Definition: DissipatingBody.cpp:422
Evolve::DissipatingBody::tidal_power
double tidal_power(unsigned zone_index, bool above, Dissipation::QuantityEntry entry=Dissipation::NO_DERIV) const
Tidal power dissipated in the given zone.
Definition: DissipatingBody.cpp:797
Evolve::DissipatingBody::angular_momentum_transfer_from_top
Eigen::Vector3d angular_momentum_transfer_from_top(unsigned zone_index, Dissipation::QuantityEntry deriv=Dissipation::NO_DERIV, bool with_respect_to_outer=false) const
Rate of angular momentum transfer (or its derivatives) to a zone due to its top boundary moving...
Definition: DissipatingBody.cpp:240
Evolve::DissipatingBody::angular_momentum_coupling
virtual Eigen::Vector3d angular_momentum_coupling(unsigned top_zone_index, Dissipation::QuantityEntry deriv=Dissipation::NO_DERIV, bool with_respect_to_top=false) const =0
Coupling torque for two neighboring zones in the coordinate system of the top zone.
Evolve::DissipatingZone::add_to_evolution
virtual void add_to_evolution()
Appends the state defined by last configure(), to the evolution.
Definition: DissipatingZone.cpp:954
Evolve::DissipatingBody::zone
virtual const DissipatingZone & zone(unsigned zone_index) const =0
A modifiable reference to one of the body&#39;s zones.
SEMIMAJOR
SEMIMAJOR
The derivative w.r.t. the semimajor axis in AU.
Definition: DissipationQuantities.h:99
MOMENT_OF_INERTIA
MOMENT_OF_INERTIA
Definition: DissipationQuantities.h:99
Evolve::zone_to_zone_transform
Eigen::Vector3d zone_to_zone_transform(const ZoneOrientation &from_zone, const ZoneOrientation &to_zone, const Eigen::Vector3d &vector, Dissipation::QuantityEntry deriv, bool with_respect_to_from)
Transforms a vector betwen the coordinates systems of two zones.
Definition: ZoneOrientation.cpp:6
Evolve::DissipatingBody::surface_lock_frequency
double surface_lock_frequency() const
Angular velocity of the surface zone when locked (assumed constant).
Definition: DissipatingBody.h:545
Evolve
Orientations of zones of bodies in a binary system.
Definition: ConstSolutionIterator.cpp:3
Evolve::DissipatingZone::tidal_torque_y
double tidal_torque_y(bool above, Dissipation::QuantityEntry entry=Dissipation::NO_DERIV) const
The dimensionless torque along y.
Definition: DissipatingZone.h:686
NUM_DERIVATIVES
NUM_DERIVATIVES
The total number of derivatives supported.
Definition: DissipationQuantities.h:99
Evolve::DissipatingZone::moment_of_inertia
virtual double moment_of_inertia(int deriv_order=0) const =0
Moment of inertia of the zone or its age derivative at the age of last configure() call...
Evolve::DissipatingBody::mass
double mass() const
The mass of the body (constant with age).
Definition: DissipatingBody.h:534
Evolve::DissipatingBody::calculate_orbit_rate_corrections
void calculate_orbit_rate_corrections()
Calculates the correction the orbital energy gain and torque due to switching locked zones from fully...
Definition: DissipatingBody.cpp:164
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::DissipatingBody::stopping_conditions
virtual CombinedStoppingCondition * stopping_conditions(BinarySystem &system, bool primary)
Conditions detecting the next possible discontinuities in the evolution due to this body...
Definition: DissipatingBody.cpp:1046
Evolve::DissipatingBody::__dorbital_frequency_da
double __dorbital_frequency_da
The derivative of the orbital frequency w.r.t. semimajor axis.
Definition: DissipatingBody.h:43
Evolve::DissipatingBody::radius
double radius(int deriv_order=0) const
The current radius or its derivative with age of the body.
Definition: DissipatingBody.h:527
Evolve::DissipatingBody::nontidal_torque
Eigen::Vector3d nontidal_torque(unsigned zone_index, Dissipation::QuantityEntry deriv=Dissipation::NO_DERIV, int deriv_zone=0) const
External torque acting on a single zone (last calculate_torques_power()).
Definition: DissipatingBody.cpp:742
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
Evolve::DissipatingZone::reset_evolution
virtual void reset_evolution()
Discards all evolution.
Definition: DissipatingZone.cpp:1004
Evolve::DissipatingBody::__orbital_frequency
double __orbital_frequency
The mean angular velocity with which the orbit is traversed.
Definition: DissipatingBody.h:43
Evolve::DissipatingBody::normalize_torques
double normalize_torques(double companion_mass, double semimajor, double orbital_frequency)
Scales the dimensionless torques as appropriate and corrects the relevant derivatives returning the n...
Definition: DissipatingBody.cpp:6
Evolve::DissipatingBody::__tidal_torques_below
std::valarray< std::valarray< Eigen::Vector3d > > __tidal_torques_below
Tidal torques and their derivatives on each zone for spin frequency approaching potential lock from b...
Definition: DissipatingBody.h:63
Core::AstroConst::Gyr
const double Gyr
Gyr [s].
Definition: AstronomicalConstants.h:24
Evolve::DissipatingZone::outer_mass
virtual double outer_mass(int deriv_order=0) const =0
Mass coordinate of the zone&#39;s outer ouboundary or its derivative.
Evolve::DissipatingBody::correct_orbit_torque
void correct_orbit_torque(std::valarray< Eigen::VectorXd > &above_lock_fractions)
Corrects __orbit_torque for zone locks.
Definition: DissipatingBody.cpp:523
ECCENTRICITY
ECCENTRICITY
The derivative w.r.t. the eccentricity.
Definition: DissipationQuantities.h:88
END_DIMENSIONLESS_DERIV
END_DIMENSIONLESS_DERIV
Definition: DissipationQuantities.h:88
Evolve::DissipatingBody::change_e_order
virtual void change_e_order(unsigned new_e_order, BinarySystem &system, bool primary)
Change the eccentricity expansion order for all dissipative zones.
Definition: DissipatingBody.cpp:1059
Core::AstroConst::day
const double day
day [s]
Definition: AstronomicalConstants.h:25
Core::AstroConst::solar_mass
const double solar_mass
Mass of the sun [kg].
Definition: AstronomicalConstants.h:21
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
Evolve::DissipatingBody::__orbit_torque_correction
std::vector< Eigen::Vector3d > __orbit_torque_correction
Corrections to __orbit_torque_below (undifferentiated) if single zones switch to above.
Definition: DissipatingBody.h:101
Evolve::DissipatingBody::angular_momentum_loss
virtual double angular_momentum_loss(Dissipation::QuantityEntry deriv=Dissipation::NO_DERIV) const =0
Rate of angular momentum loss by the top zone of the body and its derivatives.
Evolve::DissipatingBody::angular_momentum_transfer
void angular_momentum_transfer(const DissipatingZone &outer_zone, const DissipatingZone &inner_zone, Eigen::Vector3d &outer_angmom_gain, Eigen::Vector3d &inner_angmom_gain, Dissipation::QuantityEntry deriv=Dissipation::NO_DERIV, bool with_respect_to_outer=false) const
Angular momentum transfer between two neighboring zones due to moving zone boundaries.
Definition: DissipatingBody.cpp:187
Evolve::DissipatingZone::stopping_conditions
virtual CombinedStoppingCondition * stopping_conditions(BinarySystem &system, bool primary, unsigned zone_index)
Conditions detecting the next possible discontinuities in the evolution due to this zone...
Definition: DissipatingZone.cpp:1022
Evolve::DissipatingZone::locked
virtual bool locked(int orbital_frequency_multiplier, int spin_frequency_multiplier) const
Should return true iff the given term is presently locked.
Definition: DissipatingZone.h:478
NO_DERIV
NO_DERIV
The quantity itself, undifferentiated.
Definition: DissipationQuantities.h:63
DissipatingBody.h
Declares the DissipatingBody class.
ORBITAL_FREQUENCY
ORBITAL_FREQUENCY
The derivative w.r.t. the orbital frequency.
Definition: DissipationQuantities.h:63
Evolve::DissipatingZone::tidal_torque_z
double tidal_torque_z(bool above, Dissipation::QuantityEntry entry=Dissipation::NO_DERIV) const
The dimensionless tidal torque along z.
Definition: DissipatingZone.h:731
Evolve::DissipatingBody::__orbit_entries
std::vector< Dissipation::QuantityEntry > __orbit_entries
The quantities w.r.t. which derivatives of the orbit energy gain and torque are pre-calculated or err...
Definition: DissipatingBody.h:80
SPIN_ANGMOM
SPIN_ANGMOM
The derivative w.r.t. the spin angular momentum in .
Definition: DissipationQuantities.h:99
Evolve::DissipatingBody::tidal_orbit_power
double tidal_orbit_power(Dissipation::QuantityEntry entry=Dissipation::NO_DERIV, unsigned deriv_zone_index=0, const Eigen::VectorXd &above_lock_fraction_deriv=Eigen::VectorXd()) const
Rate of increase of the orbital energy due to tides in this body (last calculate_torques_power()).
Definition: DissipatingBody.cpp:849
Evolve::DissipatingZone::rewind_evolution
virtual void rewind_evolution(unsigned nsteps)
Discards the last steps from the evolution.
Definition: DissipatingZone.cpp:1012
Evolve::DissipatingBody::tidal_torque
const Eigen::Vector3d & tidal_torque(unsigned zone_index, bool above, Dissipation::QuantityEntry entry=Dissipation::NO_DERIV) const
Tidal torque acting on the given zone (last calculate_torques_power()).
Definition: DissipatingBody.h:376
Evolve::DissipatingBody::angular_momentum_transfer_to_zone
Eigen::Vector3d angular_momentum_transfer_to_zone(unsigned zone_index, Dissipation::QuantityEntry deriv=Dissipation::NO_DERIV, int deriv_zone=0) const
Rate of angular momentum transfer (or its derivatives) to a zone due to moving zone boundaries...
Definition: DissipatingBody.cpp:372
Evolve::DissipatingBody::collect_orbit_rates
void collect_orbit_rates(double orbital_frequency, double torque_norm)
Calcuates the total energy and angular momentum loss rates for the orbit due to this body&#39;s tidal dis...
Definition: DissipatingBody.cpp:79
Evolve::DissipatingBody::__tidal_torques_above
std::valarray< std::valarray< Eigen::Vector3d > > __tidal_torques_above
Tidal torques and their derivatives on each zone for spin frequency approaching potential lock from a...
Definition: DissipatingBody.h:63
Evolve::DissipatingZone::change_e_order
virtual void change_e_order(unsigned new_e_order, BinarySystem &system, bool primary, unsigned zone_index)
Changes the order of the eccentricity expansion performed.
Definition: DissipatingZone.cpp:902
Evolve::CombinedStoppingCondition
A class combining the the outputs of multiple stopping conditions.
Definition: CombinedStoppingCondition.h:21
Evolve::DissipatingBody::correct_orbit_power
void correct_orbit_power(Eigen::VectorXd &above_lock_fractions_age_deriv, Eigen::VectorXd &above_lock_fractions_semimajor_deriv, Eigen::VectorXd &above_lock_fractions_eccentricity_deriv, Eigen::VectorXd &above_lock_fractions_radius_deriv)
Corrects __orbit_power for zone locks.
Definition: DissipatingBody.cpp:438
Evolve::DissipatingBody::__orbit_torque
std::vector< Eigen::Vector3d > __orbit_torque
Torque on the orbit (and its derivatives and error) due to tides on this body.
Definition: DissipatingBody.h:101
Evolve::DissipatingZone::dissipative
virtual bool dissipative() const =0
Should return true iff the zone has some non-zero dissipation.
Evolve::DissipatingBody::__above_lock_fractions
std::valarray< Eigen::VectorXd > __above_lock_fractions
The fractional contribution of the above the lock rates for locked zones and their derivatives...
Definition: DissipatingBody.h:115
Evolve::DissipatingBody::__num_locked_zones
unsigned __num_locked_zones
The number of zones currently in a spin-orbit lock.
Definition: DissipatingBody.h:111
Evolve::DissipatingZone::tidal_torque_x
double tidal_torque_x(bool above, Dissipation::QuantityEntry entry=Dissipation::NO_DERIV) const
The dimensionless tidal torque along x.
Definition: DissipatingZone.h:639
Evolve::BinarySystem
Describes a system of two bodies orbiting each other.
Definition: BinarySystem.h:56
Evolve::DissipatingBody::__orbit_power
std::valarray< double > __orbit_power
Total tidal power (and derivatives and error) gained by the orbit from this body. ...
Definition: DissipatingBody.h:88
Evolve::DissipatingBody::__orbit_power_correction
std::valarray< double > __orbit_power_correction
Corrections to __orbit_power_below (undifferentiated) if single zones switch to above.
Definition: DissipatingBody.h:88
Evolve::DissipatingBody::tidal_orbit_torque
Eigen::Vector3d tidal_orbit_torque(Dissipation::QuantityEntry deriv=Dissipation::NO_DERIV, unsigned deriv_zone_index=0, const Eigen::VectorXd &above_lock_fraction_deriv=Eigen::VectorXd()) const
The torque on the orbit due to tidal dissipation in the body.
Definition: DissipatingBody.cpp:906
Evolve::DissipatingZone::outer_radius
virtual double outer_radius(int deriv_order=0) const =0
Outer radius of the zone or its derivative (per last.
Evolve::DissipatingBody::spin_frequency
double spin_frequency() const
The surface spin freuqency of the body.
Definition: DissipatingBody.h:538
Evolve::DissipatingZone::spin_frequency
double spin_frequency() const
The spin frequency of the given zone.
Definition: DissipatingZone.h:584
Core::AstroConst::G
const double G
Gravitational constant in SI.
Definition: AstronomicalConstants.h:26
Evolve::DissipatingBody::number_zones
virtual unsigned number_zones() const =0
The number of zones the body consists of.