Planetary Orbital Evolution due to Tides
Orbital evolution of two objects experiencing tides
BinarySystem.cpp
Go to the documentation of this file.
1 
8 #define BUILDING_LIBRARY
9 #include "BinarySystem.h"
10 
11 namespace Evolve {
12 
13  void BinarySystem::find_locked_zones()
14  {
15  __locked_zones.clear();
16  for(short body_ind = 0; body_ind < 2; ++body_ind) {
17  DissipatingBody &body = (body_ind == 0 ? __body1 : __body2);
18  for(
19  unsigned zone_ind = 0;
20  zone_ind < body.number_zones();
21  ++zone_ind
22  )
23  if(body.zone(zone_ind).locked()) {
24  body.zone(zone_ind).locked_zone_index() =
25  __locked_zones.size();
26  __locked_zones.push_back(zone_ind);
27  }
28  }
29  }
30 
31  int BinarySystem::locked_surface_differential_equations(
32  double *evolution_rates,
33  bool expansion_error
34  ) const
35  {
36  if(!expansion_error) {
37  DissipatingZone &locked_zone = __body1.zone(0);
38  locked_zone.set_evolution_rates(
39  locked_zone.moment_of_inertia(1) * locked_zone.spin_frequency(),
40  0.0,
41  0.0
42  );
43  }
44  for(
45  unsigned zone_index = 1;
46  zone_index < __body1.number_zones();
47  ++zone_index
48  ) {
49  evolution_rates[zone_index - 1] = (
50  expansion_error
51  ? 0.0
52  : __body1.nontidal_torque(zone_index)[2]
53  );
54  if(!expansion_error)
55  __body1.zone(zone_index).set_evolution_rates(
56  evolution_rates[zone_index - 1],
57  0.0,
58  0.0
59  );
60 
61  assert(__body1.nontidal_torque(zone_index)[0] == 0);
62  assert(__body1.nontidal_torque(zone_index)[1] == 0);
63 
64  }
65  return 0;
66  }
67 
68  void BinarySystem::locked_surface_jacobian(double *param_derivs,
69  double *age_derivs) const
70  {
71  unsigned num_param = __body1.number_zones() - 1;
72  double dR_dt = __body1.radius(1),
73  dIabove_dt = __body1.zone(0).moment_of_inertia(1),
74  dI_dt = __body1.zone(1).moment_of_inertia(1);
75  for(unsigned row = 0; row < num_param; ++row) {
76  double dIbelow_dt = 0;
77  age_derivs[row] =
78  __body1.nontidal_torque(row + 1, Dissipation::AGE)[2]
79  +
80  __body1.nontidal_torque(row + 1, Dissipation::RADIUS)[2] * dR_dt
81  +
82  __body1.nontidal_torque(row + 1,
83  Dissipation::MOMENT_OF_INERTIA,
84  0)[2] * dI_dt
85  +
86  __body1.nontidal_torque(row + 1,
87  Dissipation::MOMENT_OF_INERTIA,
88  -1)[2] * dIabove_dt;
89  if(row < num_param - 1) {
90  dIbelow_dt = __body1.zone(row + 2).moment_of_inertia(1);
91  age_derivs[row] += __body1.nontidal_torque(
92  row + 1,
93  Dissipation::MOMENT_OF_INERTIA,
94  1)[2] * dIbelow_dt;
95  }
96  dIabove_dt = dI_dt;
97  dI_dt = dIbelow_dt;
98 
99  assert(__body1.nontidal_torque(row + 1, Dissipation::AGE)[0] == 0);
100  assert(__body1.nontidal_torque(row + 1, Dissipation::AGE)[1] == 0);
101  assert(__body1.nontidal_torque(row + 1, Dissipation::RADIUS)[0]
102  ==
103  0);
104  assert(__body1.nontidal_torque(row + 1, Dissipation::RADIUS)[1]
105  ==
106  0);
107  assert(__body1.nontidal_torque(row + 1,
108  Dissipation::MOMENT_OF_INERTIA,
109  -1)[0]
110  ==
111  0);
112  assert(__body1.nontidal_torque(row + 1,
113  Dissipation::MOMENT_OF_INERTIA,
114  -1)[1]
115  ==
116  0);
117  assert(__body1.nontidal_torque(row + 1,
118  Dissipation::MOMENT_OF_INERTIA,
119  0)[0]
120  ==
121  0);
122  assert(__body1.nontidal_torque(row + 1,
123  Dissipation::MOMENT_OF_INERTIA,
124  0)[1]
125  ==
126  0);
127  assert(__body1.nontidal_torque(row + 1,
128  Dissipation::MOMENT_OF_INERTIA,
129  1)[0]
130  ==
131  0);
132  assert(__body1.nontidal_torque(row + 1,
133  Dissipation::MOMENT_OF_INERTIA,
134  1)[1]
135  ==
136  0);
137 
138  for(unsigned col = 0; col < num_param; ++col) {
139  double &dest = *(param_derivs + row * num_param + col);
140  if(std::abs(static_cast<int>(row) - static_cast<int>(col)) < 1)
141  dest = 0;
142  else {
143  dest = __body1.nontidal_torque(row + 1,
144  Dissipation::SPIN_ANGMOM,
145  col - row)[2];
146 
147  assert(__body1.nontidal_torque(row + 1,
148  Dissipation::SPIN_ANGMOM,
149  col - row)[0]
150  ==
151  0);
152  assert(__body1.nontidal_torque(row + 1,
153  Dissipation::SPIN_ANGMOM,
154  col - row)[1]
155  ==
156  0);
157 
158  }
159  }
160  }
161  }
162 
163  int BinarySystem::single_body_differential_equations(
164  double *evolution_rates,
165  bool expansion_error
166  ) const
167  {
168  unsigned nzones = __body1.number_zones();
169  double ref_angmom = __body1.zone(0).angular_momentum(),
170  *inclination_evol = evolution_rates,
171  *periapsis_evol = evolution_rates + (nzones - 1),
172  *angmom_evol = periapsis_evol+(nzones - 1);
173 
174  if(expansion_error) {
175  angmom_evol[0] = 0.0;
176  for(
177  unsigned zone_index = 0;
178  zone_index < nzones - 1;
179  ++zone_index
180  )
181  inclination_evol[zone_index] = (
182  periapsis_evol[zone_index] = (
183  angmom_evol[zone_index + 1] = 0.0
184  )
185  );
186  return 0;
187  }
188 
189  Eigen::Vector3d reference_torque;
190  for(unsigned zone_index = 0; zone_index < nzones; ++zone_index) {
191  Eigen::Vector3d torque = __body1.nontidal_torque(zone_index);
192  angmom_evol[zone_index] = torque[2];
193 
194  DissipatingZone &zone = __body1.zone(zone_index);
195  if(zone_index) {
196  zone.set_reference_zone_angmom(ref_angmom);
197  inclination_evol[zone_index - 1] = zone.inclination_evolution(
198  zone_to_zone_transform(__body1.zone(0),
199  zone,
200  reference_torque),
201  torque
202  );
203  periapsis_evol[zone_index - 1] = zone.periapsis_evolution(
204  zone_to_zone_transform(__body1.zone(0),
205  zone,
206  reference_torque),
207  torque
208  );
209  } else reference_torque = torque;
210 
211  if(!expansion_error)
212  zone.set_evolution_rates(
213  angmom_evol[zone_index],
214  (zone_index ? inclination_evol[zone_index - 1] : 0.0),
215  (zone_index ? periapsis_evol[zone_index -1] : 0.0)
216  );
217  }
218 #ifdef VERBOSE_DEBUG
219  std::cerr << "rates: ";
220  for(unsigned i = 0; i < 3 * nzones - 2; ++i) {
221  if(i) std::cerr << ", ";
222  std::cerr << evolution_rates[i];
223  }
224  std::cerr << std::endl;
225 #endif
226  return 0;
227  }
228 
229  void BinarySystem::fill_single_body_jacobian(
230  double *inclination_param_derivs,
231  double *periapsis_param_derivs,
232  double *angmom_param_derivs,
233  double *inclination_age_derivs,
234  double *periapsis_age_derivs,
235  double *angmom_age_derivs
236  ) const
237  {
238  unsigned nzones=__body1.number_zones(), nparams = 3 * nzones - 2;
239  Eigen::Vector3d ref_torque = __body1.nontidal_torque(0),
240  dref_torque_dincl = __body1.nontidal_torque(
241  0,
242  Dissipation::INCLINATION,
243  1
244  ),
245  dref_torque_dperi=__body1.nontidal_torque(
246  0,
247  Dissipation::PERIAPSIS,
248  1
249  ),
250  dref_torque_dangmom0=__body1.nontidal_torque(
251  0,
252  Dissipation::SPIN_ANGMOM,
253  0
254  ),
255  dref_torque_dangmom1=__body1.nontidal_torque(
256  0,
257  Dissipation::SPIN_ANGMOM,
258  1
259  ),
260  dref_torque_dage = __body1.nontidal_torque(
261  0,
262  Dissipation::AGE
263  ),
264  zero3d(0, 0, 0);
265  for(
266  unsigned deriv_zone_ind = 0;
267  deriv_zone_ind < nzones;
268  ++deriv_zone_ind
269  ) {
270  angmom_param_derivs[deriv_zone_ind - 1] =
271  (deriv_zone_ind == 1 ? dref_torque_dincl[2] : 0);
272  angmom_param_derivs[nzones+deriv_zone_ind - 2] =
273  (deriv_zone_ind == 1 ? dref_torque_dperi[2] : 0);
274  double &angmom_angmom_deriv =
275  angmom_param_derivs[2 * nzones+deriv_zone_ind - 2];
276  if(deriv_zone_ind==0) angmom_angmom_deriv = dref_torque_dangmom0[2];
277  else if(deriv_zone_ind == 1)
278  angmom_angmom_deriv = dref_torque_dangmom1[2];
279  else angmom_angmom_deriv = 0;
280  }
281  angmom_age_derivs[0] = dref_torque_dage[2];
282  DissipatingZone &ref_zone = __body1.zone(0);
283  for(unsigned zone_ind = 1; zone_ind < nzones; ++zone_ind) {
284  DissipatingZone &zone = __body1.zone(zone_ind);
285  Eigen::Vector3d
286  zone_ref_torque = zone_to_zone_transform(ref_zone,
287  zone,
288  ref_torque),
289  zone_torque = __body1.nontidal_torque(zone_ind),
290  ref_torque_age_deriv = zone_to_zone_transform(ref_zone,
291  zone,
292  dref_torque_dage),
293  zone_torque_age_deriv=__body1.nontidal_torque(zone_ind,
294  Dissipation::AGE,
295  0);
296  inclination_age_derivs[zone_ind - 1] =
297  zone.inclination_evolution(zone_ref_torque,
298  zone_torque,
299  Dissipation::AGE,
300  ref_torque_age_deriv,
301  zone_torque_age_deriv);
302  periapsis_age_derivs[zone_ind - 1] =
303  zone.periapsis_evolution(zone_ref_torque,
304  zone_torque,
305  Dissipation::AGE,
306  ref_torque_age_deriv,
307  zone_torque_age_deriv);
308  angmom_age_derivs[zone_ind] = zone_torque_age_deriv[2];
309  for(
310  unsigned quantity_ind = 0;
311  quantity_ind < nparams;
312  ++quantity_ind
313  ) {
314  unsigned offset = (zone_ind - 1) * nparams + quantity_ind;
315  double &inclination_dest = inclination_param_derivs[offset],
316  &periapsis_dest = periapsis_param_derivs[offset],
317  &angmom_dest = angmom_param_derivs[offset+nparams];
318  unsigned quantity_zone = quantity_ind;
319  Dissipation::QuantityEntry with_respect_to;
320  Eigen::Vector3d ref_torque_deriv(0, 0, 0);
321  if(quantity_ind > 2 * nzones - 2) {
322  quantity_zone -= 2 * nzones - 2;
323  with_respect_to = Dissipation::SPIN_ANGMOM;
324  if(quantity_zone == 0)
325  ref_torque_deriv = dref_torque_dangmom0;
326  else if(quantity_zone == 1)
327  ref_torque_deriv = dref_torque_dangmom1;
328  if(quantity_zone <= 1)
329  ref_torque_deriv = zone_to_zone_transform(
330  ref_zone,
331  zone,
332  ref_torque_deriv
333  );
334  } else {
335  if(quantity_ind >= nzones - 1) {
336  quantity_zone -= nzones - 2;
337  with_respect_to = Dissipation::PERIAPSIS;
338  if(quantity_zone == 1)
339  ref_torque_deriv = dref_torque_dperi;
340  } else {
341  quantity_zone += 1;
342  with_respect_to = Dissipation::INCLINATION;
343  if(quantity_zone == 1)
344  ref_torque_deriv = dref_torque_dincl;
345  }
346  if(quantity_zone == 1)
347  ref_torque_deriv = zone_to_zone_transform(
348  ref_zone,
349  zone,
350  ref_torque_deriv
351  );
352  ref_torque_deriv += zone_to_zone_transform(ref_zone,
353  zone,
354  ref_torque,
355  with_respect_to);
356  }
357  Eigen::Vector3d zone_torque_deriv;
358  if(quantity_zone == zone_ind) {
359  zone_torque_deriv = __body1.nontidal_torque(zone_ind,
360  with_respect_to,
361  0);
362  inclination_dest=zone.inclination_evolution(
363  zone_ref_torque,
364  zone_torque,
365  with_respect_to,
366  ref_torque_deriv,
367  zone_torque_deriv
368  );
369  periapsis_dest=zone.periapsis_evolution(zone_ref_torque,
370  zone_torque,
371  with_respect_to,
372  ref_torque_deriv,
373  zone_torque_deriv);
374  } else {
375  if(std::abs(static_cast<int>(quantity_zone-zone_ind)) == 1) {
376  zone_torque_deriv = __body1.nontidal_torque(
377  zone_ind,
378  with_respect_to,
379  quantity_zone - zone_ind
380  );
381  } else zone_torque_deriv = zero3d;
382  inclination_dest = zone.inclination_evolution(
383  ref_torque_deriv,
384  zone_torque_deriv
385  );
386  periapsis_dest = zone.periapsis_evolution(ref_torque_deriv,
387  zone_torque_deriv);
388  }
389  angmom_dest = zone_torque_deriv[2];
390  }
391  }
392  }
393 
394  void BinarySystem::single_body_jacobian(double *param_derivs,
395  double *age_derivs) const
396  {
397  unsigned nzones = __body1.number_zones(),
398  nparams=3*nzones-2;
399  fill_single_body_jacobian(param_derivs,
400  param_derivs + (nzones - 1) * nparams,
401  param_derivs + 2 * (nzones - 1) * nparams,
402  age_derivs,
403  age_derivs + (nzones - 1),
404  age_derivs + 2 * (nzones - 1));
405  }
406 
407  double BinarySystem::semimajor_evolution(
408  double orbit_power,
409  double orbit_power_deriv
410  ) const
411  {
412  if(std::isnan(orbit_power_deriv))
413  return (orbit_power == 0
414  ? 0
415  : -__semimajor * orbit_power / __orbital_energy);
416 
417  else if(orbit_power == 0 && orbit_power_deriv == 0)
418  return 0;
419 
420  return (
421  -(2.0 * orbit_power
422  +
423  __semimajor * orbit_power_deriv
424  )
425  /
426  __orbital_energy
427  );
428  }
429 
430  double BinarySystem::eccentricity_evolution(double orbit_power,
431  double orbit_angmom_gain,
432  double orbit_power_deriv,
433  double orbit_angmom_gain_deriv,
434  bool semimajor_deriv) const
435  {
436  if(__eccentricity <= 1e-8) return 0;
437  double e2 = std::pow(__eccentricity, 2),
438  factor = -(1.0 - e2) / (2.0 * __eccentricity);
439 
440  if(std::isnan(orbit_power_deriv))
441  return factor * (orbit_power / __orbital_energy
442  +
443  2.0 * orbit_angmom_gain / __orbital_angmom);
444  else if(semimajor_deriv)
445  return (
446  factor
447  *
448  (
449  (
450  orbit_power_deriv
451  +
452  orbit_power / __semimajor
453  )
454  /
455  __orbital_energy
456  +
457  (
458  2.0 * orbit_angmom_gain_deriv
459  -
460  orbit_angmom_gain / __semimajor
461  )
462  /
463  __orbital_angmom
464  )
465  );
466  else return (
467  factor
468  *
469  (
470  orbit_power_deriv / __orbital_energy
471  +
472  2.0 * orbit_angmom_gain_deriv / __orbital_angmom
473  )
474  -
475  2.0 * orbit_angmom_gain / __orbital_angmom
476  -
477  (1.0 + e2) / e2 * (
478  orbit_power / __orbital_energy
479  +
480  2.0 * orbit_angmom_gain / __orbital_angmom
481  )
482  );
483  }
484 
485  double BinarySystem::eccentricity_evolution_expansion_error() const
486  {
487  if(__eccentricity <= 1e-8) return 0;
488  double factor = ((1.0 - std::pow(__eccentricity, 2))
489  /
490  (2.0 * __eccentricity));
491 
492  return factor * (
493  std::abs(__orbit_power_expansion_error / __orbital_energy)
494  +
495  std::abs(2.0 * __orbit_angmom_gain_expansion_error
496  /
497  __orbital_angmom)
498  );
499  }
500 
501  void BinarySystem::above_lock_problem_deriv_correction(
502  Dissipation::QuantityEntry entry,
503  bool body1_deriv,
504  Eigen::MatrixXd &matrix,
505  Eigen::VectorXd &rhs
506  ) const
507  {
508  unsigned deriv_zone_index = (body1_deriv
509  ? 0
510  : __body1.number_locked_zones());
511  DissipatingBody &deriv_body = (body1_deriv ? __body1 : __body2);
512  DissipatingZone &deriv_zone = deriv_body.zone(0);
513 
514  if(
515  entry == Dissipation::ORBITAL_FREQUENCY
516  ||
517  entry == Dissipation::SEMIMAJOR
518  ) {
519  double coef = (entry == Dissipation::SEMIMAJOR
520  ? 1
521  : Core::orbital_angular_velocity(__body1.mass(),
522  __body2.mass(),
523  __semimajor,
524  true));
525  bool done = false;
526  unsigned locked_ind = 0;
527  for(DissipatingBody *body = &__body1; !done; body = &__body2) {
528  for(
529  unsigned zone_ind = 0;
530  zone_ind < body->number_zones();
531  ++zone_ind
532  ) {
533  rhs.array() += (1.5
534  *
535  (
536  __body1.tidal_orbit_power()
537  +
538  __body2.tidal_orbit_power()
539  )
540  /
541  __orbital_energy
542  /
543  __semimajor
544  *
545  coef);
546  if(body->zone(zone_ind).locked()) {
547  matrix.col(locked_ind).array() += (
548  1.5
549  *
550  (
551  body->tidal_power(zone_ind, true)
552  -
553  body->tidal_power(zone_ind, false)
554  )
555  /
556  __orbital_energy
557  /
558  __semimajor * coef
559  );
560 
561  ++locked_ind;
562  }
563  }
564  done = (body == &__body2);
565  }
566  } else if(deriv_zone.locked()) {
567  if(
568  entry == Dissipation::SPIN_FREQUENCY
569  ||
570  entry == Dissipation::SPIN_ANGMOM
571  ) {
572  double coef = (
573  entry == Dissipation::SPIN_FREQUENCY
574  ? deriv_zone.moment_of_inertia()
575  : 1
576  ) / std::pow(deriv_zone.angular_momentum(), 2);
577  matrix(deriv_zone_index, deriv_zone_index) -=(
578  deriv_body.tidal_torque(0, true)[2]
579  -
580  deriv_body.tidal_torque(0, false)[2]
581  ) * coef;
582  rhs(deriv_zone_index) -= (
583  deriv_body.tidal_torque(0, false)[2]
584  *
585  coef
586  );
587  } else if(entry == Dissipation::MOMENT_OF_INERTIA) {
588  rhs(deriv_zone_index) -= (
589  deriv_zone.moment_of_inertia(1)
590  /
591  std::pow(deriv_zone.moment_of_inertia(), 2)
592  );
593  }
594  }
595  }
596 
597  void BinarySystem::calculate_above_lock_fractions(
598  Eigen::VectorXd &fractions,
599  Dissipation::QuantityEntry entry,
600  bool body1_deriv
601  )
602  {
603  unsigned num_locked_zones = __locked_zones.size();
604 
605  assert(num_locked_zones == (__body1.number_locked_zones()
606  +
607  __body2.number_locked_zones()));
608 
609  if(num_locked_zones == 0) {
610  fractions.resize(0);
611  return;
612  }
613  std::valarray<double> nontidal_torque(num_locked_zones),
614  tidal_torque_z_above(num_locked_zones),
615  tidal_torque_z_below(num_locked_zones),
616  tidal_power_difference(num_locked_zones);
617  unsigned locked_zone_ind = 0;
618  for(
619  std::list<unsigned>::const_iterator zi = __locked_zones.begin();
620  zi != __locked_zones.end();
621  ++zi
622  ) {
623  DissipatingBody &body = (
624  locked_zone_ind < __body1.number_locked_zones()
625  ? __body1
626  : __body2
627  );
628  nontidal_torque[locked_zone_ind] =
629  body.nontidal_torque(*zi, entry)[2];
630  tidal_torque_z_above[locked_zone_ind] =
631  body.tidal_torque(*zi, true, entry)[2];
632  tidal_torque_z_below[locked_zone_ind] =
633  body.tidal_torque(*zi, false, entry)[2];
634  if(!zone_specific(entry) || *zi == 0) {
635  tidal_power_difference[locked_zone_ind] = (
636  body.tidal_power(*zi, true, entry)
637  -
638  body.tidal_power(*zi, false, entry)
639  );
640  }
641  ++locked_zone_ind;
642  }
643  Eigen::MatrixXd matrix(num_locked_zones, num_locked_zones);
644  Eigen::VectorXd rhs(num_locked_zones);
645  rhs.setConstant(1.5
646  *
647  (
648  __body1.tidal_orbit_power(entry)
649  +
650  __body2.tidal_orbit_power(entry)
651  )
652  /
653  __orbital_energy);
654  unsigned i = 0;
655  for(
656  std::list<unsigned>::const_iterator
657  zi=__locked_zones.begin();
658  zi!=__locked_zones.end();
659  ++zi
660  ) {
661  if(!zone_specific(entry) || *zi == 0)
662  matrix.col(i).setConstant(1.5
663  *
664  tidal_power_difference[i]
665  /
666  __orbital_energy);
667  else matrix.col(i).setZero();
668  DissipatingZone &zone = (i < __body1.number_locked_zones()
669  ? __body1
670  : __body2).zone(*zi);
671  matrix(i, i) += (
672  (tidal_torque_z_above[i] - tidal_torque_z_below[i])
673  /
674  zone.angular_momentum()
675  );
676  if(entry == Dissipation::NO_DERIV)
677  rhs(i) += (zone.moment_of_inertia(1)
678  /
679  zone.moment_of_inertia());
680  else if(entry == Dissipation::AGE)
681  rhs(i) += (zone.moment_of_inertia(2)
682  /
683  zone.moment_of_inertia());
684  rhs(i) -= ((nontidal_torque[i] + tidal_torque_z_below[i])
685  /
686  zone.angular_momentum());
687  ++i;
688  }
689  above_lock_problem_deriv_correction(entry, body1_deriv, matrix, rhs);
690  if(entry == Dissipation::NO_DERIV) {
691  __above_lock_fractions_decomp.compute(matrix);
692  fractions = __above_lock_fractions_decomp.solve(rhs);
693  } else {
694  fractions = __above_lock_fractions_decomp.solve(
695  rhs - matrix * __above_lock_fractions[Dissipation::NO_DERIV]
696  );
697  }
698  }
699 
700  Eigen::VectorXd BinarySystem::above_lock_fractions_deriv(
701  Dissipation::QuantityEntry entry,
702  DissipatingBody &body,
703  unsigned zone_index)
704  {
705  assert(entry == Dissipation::INCLINATION
706  ||
707  entry == Dissipation::PERIAPSIS
708  ||
709  entry == Dissipation::MOMENT_OF_INERTIA
710  ||
711  entry == Dissipation::SPIN_ANGMOM);
712  assert(zone_index > 0 || &body == &__body2);
713  assert(body.number_zones() > zone_index);
714  assert(number_locked_zones()
715  ==
716  __above_lock_fractions[Dissipation::NO_DERIV].size());
717 
718  unsigned num_locked_zones = number_locked_zones();
719  if(num_locked_zones == 0) return Eigen::VectorXd();
720  DissipatingZone &deriv_zone = body.zone(zone_index);
721  unsigned locked_zone_index=(&body == &__body1
722  ? 0
723  : __body1.number_locked_zones());
724  Eigen::VectorXd rhs;
725  for(unsigned i = 0; i < zone_index; ++i)
726  if(body.zone(i).locked()) ++locked_zone_index;
727  if(deriv_zone.locked()) {
728  double above_frac =
729  __above_lock_fractions[Dissipation::NO_DERIV][locked_zone_index];
730  rhs.setConstant(
731  num_locked_zones,
732  -1.5 * (above_frac*body.tidal_power(zone_index, true, entry)
733  +
734  (1.0 - above_frac) * body.tidal_power(zone_index,
735  false,
736  entry))
737  -
738  body.nontidal_torque(zone_index, entry)[2]
739  );
740  rhs(locked_zone_index) -=
741  above_frac*body.tidal_torque(zone_index, true, entry)[2]
742  +
743  (1.0 - above_frac) * body.tidal_torque(zone_index,
744  false,
745  entry)[2];
746  if(entry == Dissipation::MOMENT_OF_INERTIA)
747  rhs(locked_zone_index) -= (
748  deriv_zone.moment_of_inertia(1)
749  /
750  std::pow(deriv_zone.moment_of_inertia(), 2)
751  );
752  else if(entry == Dissipation::SPIN_ANGMOM)
753  rhs(locked_zone_index) += (
754  above_frac*body.tidal_torque(zone_index, true)[2]
755  +
756  (1.0 - above_frac)
757  *
758  body.tidal_torque(zone_index, false)[2]
759  +
760  body.nontidal_torque(zone_index)[2]
761  ) / std::pow(deriv_zone.angular_momentum(), 2);
762  } else {
763  rhs.setConstant(num_locked_zones,
764  -1.5 * body.tidal_power(zone_index, false, entry));
765  }
766  if(zone_index > 0 && body.zone(zone_index - 1).locked())
767  rhs(body.zone(zone_index - 1).locked_zone_index()) -=
768  body.nontidal_torque(zone_index - 1, entry, 1)[2];
769  if(
770  zone_index + 1 < body.number_zones()
771  &&
772  body.zone(zone_index + 1).locked()
773  )
774  rhs(body.zone(zone_index + 1).locked_zone_index()) -=
775  body.nontidal_torque(zone_index + 1, entry, -1)[2];
776  return __above_lock_fractions_decomp.solve(rhs);
777  }
778 
779  void BinarySystem::fill_above_lock_fractions_deriv()
780  {
781  unsigned num_zones = __body1.number_zones() + __body2.number_zones();
782  DissipatingBody *body = &__body1;
783  __above_lock_fractions_inclination_deriv.resize(num_zones);
784  __above_lock_fractions_periapsis_deriv.resize(num_zones);
785  __above_lock_fractions_inertia_deriv.resize(num_zones);
786  __above_lock_fractions_angmom_deriv.resize(num_zones);
787  __above_lock_fractions_inclination_deriv[0] =
788  __above_lock_fractions[Dissipation::INCLINATION];
789  __above_lock_fractions_periapsis_deriv[0] =
790  __above_lock_fractions[Dissipation::PERIAPSIS];
791  __above_lock_fractions_inertia_deriv[0] =
792  __above_lock_fractions[Dissipation::MOMENT_OF_INERTIA];
793  __above_lock_fractions_angmom_deriv[0] =
794  __above_lock_fractions[Dissipation::SPIN_ANGMOM];
795  unsigned body_zone_ind = 0;
796  for(unsigned zone_ind = 1; zone_ind < num_zones; ++zone_ind) {
797  ++body_zone_ind;
798  if(body_zone_ind == body->number_zones()) {
799  body_zone_ind = 0;
800  body = &__body2;
801  }
802  __above_lock_fractions_inclination_deriv[zone_ind] =
803  above_lock_fractions_deriv(Dissipation::INCLINATION,
804  *body,
805  body_zone_ind);
806  __above_lock_fractions_periapsis_deriv[zone_ind] =
807  above_lock_fractions_deriv(Dissipation::PERIAPSIS,
808  *body,
809  body_zone_ind);
810  __above_lock_fractions_inertia_deriv[zone_ind] =
811  above_lock_fractions_deriv(Dissipation::MOMENT_OF_INERTIA,
812  *body,
813  body_zone_ind);
814  __above_lock_fractions_angmom_deriv[zone_ind] =
815  above_lock_fractions_deriv(Dissipation::SPIN_ANGMOM,
816  *body,
817  body_zone_ind);
818  }
819  }
820 
821  void BinarySystem::update_above_lock_fractions()
822  {
823  std::valarray<Eigen::VectorXd>
824  body2_above_lock_fractions(Dissipation::NUM_DERIVATIVES);
825  for(
826  int deriv = Dissipation::NO_DERIV;
827  deriv < Dissipation::NUM_DERIVATIVES;
828  ++deriv
829  ) {
830  calculate_above_lock_fractions(
831  __above_lock_fractions[deriv],
832  static_cast<Dissipation::QuantityEntry>(deriv)
833  );
834  if(!zone_specific(static_cast < Dissipation::QuantityEntry>(deriv)))
835  body2_above_lock_fractions[deriv] = __above_lock_fractions[deriv];
836  else calculate_above_lock_fractions(
837  body2_above_lock_fractions[deriv],
838  static_cast<Dissipation::QuantityEntry>(deriv),
839  false
840  );
841  }
842  __above_lock_fractions_body2_radius_deriv =
843  body2_above_lock_fractions[Dissipation::RADIUS];
844  __body1.set_above_lock_fractions(__above_lock_fractions);
845  __body2.set_above_lock_fractions(body2_above_lock_fractions);
846  fill_above_lock_fractions_deriv();
847  }
848 
849  void BinarySystem::fill_orbit_torque_and_power()
850  {
851  __orbit_torque = (
852  __body1.tidal_orbit_torque()
853  +
854  zone_to_zone_transform(__body2.zone(0),
855  __body1.zone(0),
856  __body2.tidal_orbit_torque())
857  );
858 
859  __orbit_torque_expansion_error.setConstant(
860  __body1.tidal_orbit_torque(Dissipation::EXPANSION_ERROR).norm()
861  +
862  __body2.tidal_orbit_torque(Dissipation::EXPANSION_ERROR).norm()
863  );
864 
865  __orbit_power = (__body1.tidal_orbit_power()
866  +
867  __body2.tidal_orbit_power());
868 
869  __orbit_power_expansion_error = (
870  __body1.tidal_orbit_power(Dissipation::EXPANSION_ERROR)
871  +
872  __body2.tidal_orbit_power(Dissipation::EXPANSION_ERROR)
873  );
874 
875  __orbit_angmom_gain = (__orbit_torque[0]
876  *
877  std::sin(__body1.zone(0).inclination())
878  +
879  __orbit_torque[2]
880  *
881  std::cos(__body1.zone(0).inclination()));
882 
883 
884  __orbit_angmom_gain_expansion_error = (
885  std::abs(__orbit_torque_expansion_error[0]
886  *
887  std::sin(__body1.zone(0).inclination()))
888  +
889  std::abs(__orbit_torque_expansion_error[2]
890  *
891  std::cos(__body1.zone(0).inclination()))
892  );
893  }
894 
895  int BinarySystem::binary_differential_equations(
896  double *differential_equations,
897  bool expansion_error
898  ) const
899  {
900  DissipatingZone &reference_zone = __body1.zone(0);
901  unsigned num_body1_zones = __body1.number_zones(),
902  num_total_zones = num_body1_zones + __body2.number_zones();
903  double *inclination_rates = differential_equations + 2,
904  *periapsis_rates = inclination_rates + num_total_zones,
905  *angmom_rates = periapsis_rates + num_total_zones - 1;
906  unsigned angmom_skipped = 0;
907  double reference_periapsis_rate = Core::NaN;
908  for(unsigned zone_ind = 0; zone_ind < num_total_zones; ++zone_ind) {
909  DissipatingBody &body = (zone_ind < num_body1_zones
910  ? __body1
911  : __body2);
912  unsigned body_zone_ind = zone_ind;
913  if(zone_ind >= num_body1_zones)
914  body_zone_ind -= num_body1_zones;
915  DissipatingZone &zone = body.zone(body_zone_ind);
916  Eigen::Vector3d zone_orbit_torque,
917  total_zone_torque;
918  if(expansion_error)
919  total_zone_torque.setZero();
920  else
921  total_zone_torque = body.nontidal_torque(body_zone_ind);
922 
923  if(expansion_error)
924  zone_orbit_torque = __orbit_torque_expansion_error;
925  else
926  zone_orbit_torque = (zone_ind
927  ? zone_to_zone_transform(reference_zone,
928  zone,
929  __orbit_torque)
930  : __orbit_torque);
931 
932  Dissipation::QuantityEntry entry = (
933  expansion_error
934  ? Dissipation::EXPANSION_ERROR
935  : Dissipation::NO_DERIV
936  );
937  if(zone.locked()) {
938  double above_frac = (__above_lock_fractions
939  [Dissipation::NO_DERIV]
940  [angmom_skipped]);
941  total_zone_torque += (
942  above_frac
943  *
944  body.tidal_torque(body_zone_ind, true, entry)
945  +
946  (1.0 - above_frac)
947  *
948  body.tidal_torque(body_zone_ind, false, entry)
949  );
950  ++angmom_skipped;
951  } else total_zone_torque += body.tidal_torque(body_zone_ind,
952  false,
953  entry);
954  inclination_rates[zone_ind] = zone.inclination_evolution(
955  zone_orbit_torque,
956  total_zone_torque,
957  entry
958  );
959  assert(!std::isnan(inclination_rates[zone_ind]));
960  if(zone_ind) {
961  periapsis_rates[zone_ind - 1] = zone.periapsis_evolution(
962  zone_orbit_torque,
963  total_zone_torque,
964  entry
965  ) - reference_periapsis_rate;
966  assert(!std::isnan(periapsis_rates[zone_ind - 1]));
967  } else {
968  reference_periapsis_rate = zone.periapsis_evolution(
969  zone_orbit_torque,
970  total_zone_torque,
971  entry
972  );
973  if(expansion_error)
974  reference_periapsis_rate = -std::abs(
975  reference_periapsis_rate
976  );
977  assert(!std::isnan(reference_periapsis_rate));
978 
979  }
980  if(!zone.locked()) {
981  angmom_rates[zone_ind - angmom_skipped] = total_zone_torque[2];
982  assert(!std::isnan(angmom_rates[zone_ind - angmom_skipped]));
983  }
984 
985  if(!expansion_error)
986  zone.set_evolution_rates(
987  total_zone_torque[2],
988  inclination_rates[zone_ind],
989  (zone_ind ? periapsis_rates[zone_ind - 1] : 0.0)
990  );
991 
992 #ifdef VERBOSE_DEBUG
993  std::cerr << "Zone " << zone_ind
994  << " torque: " << total_zone_torque
995  << " inclination rate";
996  if(expansion_error)
997  std::cerr << " error";
998  std::cerr << ": " << inclination_rates[zone_ind];
999  if(zone_ind)
1000  std::cerr << " periapsis rate: "
1001  << periapsis_rates[zone_ind - 1];
1002 
1003  std::cerr << std::endl;
1004 #endif
1005  assert(!std::isnan(total_zone_torque.sum()));
1006 
1007  }
1008  differential_equations[0] = (
1009  expansion_error
1010  ? semimajor_evolution_expansion_error()
1011  : semimajor_evolution(__orbit_power)
1012  );
1013  differential_equations[1] = (
1014  expansion_error
1015  ? eccentricity_evolution_expansion_error()
1016  : eccentricity_evolution(__orbit_power, __orbit_angmom_gain)
1017  );
1018 
1019  if(!expansion_error) {
1020  __semimajor_rate = differential_equations[0];
1021  __eccentricity_rate = differential_equations[1];
1022  }
1023 
1024  if(!angmom_skipped)
1025  differential_equations[0] *= 6.5 * std::pow(__semimajor, 5.5);
1026 
1027 #ifdef VERBOSE_DEBUG
1028  if(expansion_error)
1029  std::cerr << "rate errors: ";
1030  else
1031  std::cerr << "rates: ";
1032  for(
1033  unsigned i = 0;
1034  i < 3 * (__body1.number_zones() + __body2.number_zones()) + 1;
1035  ++i
1036  ) {
1037  if(i) std::cerr << ", ";
1038  std::cerr << differential_equations[i];
1039  }
1040  std::cerr << std::endl;
1041 #endif
1042 
1043  return 0;
1044  }
1045 
1046  template<typename VALUE_TYPE>
1047  void BinarySystem::add_body_rate_deriv(
1048  const DissipatingBody &body,
1049  VALUE_TYPE (DissipatingBody::*func)(Dissipation::QuantityEntry,
1050  unsigned,
1051  const Eigen::VectorXd &) const,
1052  std::valarray<VALUE_TYPE> &orbit_rate_deriv,
1053  unsigned offset
1054  ) const
1055  {
1056  unsigned num_zones = __body1.number_zones() + __body2.number_zones(),
1057  num_param = orbit_rate_deriv.size() - 1;
1058  orbit_rate_deriv[0] += (body.*func)(Dissipation::SEMIMAJOR,
1059  0,
1060  Eigen::VectorXd());
1061  orbit_rate_deriv[1] += (body.*func)(Dissipation::ECCENTRICITY,
1062  0,
1063  Eigen::VectorXd());
1064  orbit_rate_deriv[num_param] += (
1065  (body.*func)(Dissipation::AGE, 0, Eigen::VectorXd())
1066  +
1067  (body.*func)(Dissipation::RADIUS, 0, Eigen::VectorXd())
1068  *
1069  body.radius(1)
1070  );
1071  unsigned locked_zone_count = (offset ? __body1.number_locked_zones() : 0);
1072  for(unsigned zone_ind = 0; zone_ind < body.number_zones(); ++zone_ind) {
1073  orbit_rate_deriv[zone_ind+2+offset] =
1074  (body.*func)(
1075  Dissipation::INCLINATION,
1076  zone_ind,
1077  __above_lock_fractions_inclination_deriv[zone_ind + offset]
1078  );
1079  if(zone_ind+offset > 0)
1080  orbit_rate_deriv[zone_ind + 1 + num_zones + offset] =
1081  (body.*func)(
1082  Dissipation::PERIAPSIS,
1083  zone_ind,
1084  __above_lock_fractions_periapsis_deriv[zone_ind+offset]
1085  );
1086  const DissipatingZone &zone = body.zone(zone_ind);
1087  if(zone.locked()) ++locked_zone_count;
1088  else orbit_rate_deriv[
1089  zone_ind + 1 + 2 * num_zones - locked_zone_count
1090  ] = (body.*func)(
1091  Dissipation::SPIN_ANGMOM,
1092  zone_ind,
1093  __above_lock_fractions_angmom_deriv[zone_ind + offset]
1094  );
1095  orbit_rate_deriv[num_param] += (
1096  (body.*func)(
1097  Dissipation::MOMENT_OF_INERTIA,
1098  zone_ind,
1099  __above_lock_fractions_inertia_deriv[zone_ind + offset]
1100  )
1101  *
1102  zone.moment_of_inertia(1)
1103  );
1104  }
1105  }
1106 
1107  void BinarySystem::fill_orbit_power_deriv(
1108  std::valarray<double> &orbit_power_deriv) const
1109  {
1110  orbit_power_deriv[0] = 0;
1111  orbit_power_deriv[1] = 0;
1112  orbit_power_deriv[orbit_power_deriv.size() - 1] = 0;
1113  add_body_rate_deriv(__body1, &DissipatingBody::tidal_orbit_power,
1114  orbit_power_deriv,
1115  0);
1116  add_body_rate_deriv(__body2,
1117  &DissipatingBody::tidal_orbit_power,
1118  orbit_power_deriv,
1119  __body1.number_zones());
1120  }
1121 
1122  void BinarySystem::fill_orbit_angmom_gain_deriv(
1123  std::valarray<double> &orbit_angmom_gain_deriv
1124  ) const
1125  {
1126  std::valarray<Eigen::Vector3d>
1127  body1_orbit_torque_deriv(orbit_angmom_gain_deriv.size()),
1128  body2_orbit_torque_deriv(orbit_angmom_gain_deriv.size());
1129  body1_orbit_torque_deriv[0] =
1130  body1_orbit_torque_deriv[1] =
1131  body2_orbit_torque_deriv[0] =
1132  body2_orbit_torque_deriv[1] = Eigen::Vector3d(0, 0, 0);
1133  add_body_rate_deriv(__body1,
1134  &DissipatingBody::tidal_orbit_torque,
1135  body1_orbit_torque_deriv,
1136  0);
1137  add_body_rate_deriv(__body2, &DissipatingBody::tidal_orbit_torque,
1138  body2_orbit_torque_deriv,
1139  __body1.number_zones());
1140  double body1_sin_inc = std::sin(__body1.zone(0).inclination()),
1141  body1_cos_inc = std::sin(__body1.zone(0).inclination()),
1142  body2_sin_inc = std::sin(__body2.zone(0).inclination()),
1143  body2_cos_inc = std::sin(__body2.zone(0).inclination());
1144  for(unsigned i = 0; i < orbit_angmom_gain_deriv.size(); ++i)
1145  orbit_angmom_gain_deriv[i] = (
1146  body1_sin_inc * body1_orbit_torque_deriv[i][0]
1147  +
1148  body1_cos_inc * body1_orbit_torque_deriv[i][2]
1149  +
1150  body2_sin_inc * body2_orbit_torque_deriv[i][0]
1151  +
1152  body2_cos_inc * body2_orbit_torque_deriv[i][2]
1153  );
1154  Eigen::Vector3d body1_orbit_torque = __body1.tidal_orbit_torque(),
1155  body2_orbit_torque = __body2.tidal_orbit_torque();
1156  orbit_angmom_gain_deriv[2] +=(body1_cos_inc * body1_orbit_torque[0]
1157  -
1158  body1_sin_inc * body1_orbit_torque[2]
1159  +
1160  body2_cos_inc * body2_orbit_torque[0]
1161  -
1162  body2_sin_inc * body2_orbit_torque[2]);
1163  }
1164 
1165  void BinarySystem::semimajor_jacobian(
1166  const std::valarray<double> &orbit_power_deriv,
1167  bool a6p5,
1168  double *param_derivs,
1169  double &age_deriv
1170  ) const
1171  {
1172  param_derivs[0] = semimajor_evolution(__orbit_power,
1173  orbit_power_deriv[0]);
1174  if(a6p5) param_derivs[0] +=
1175  5.5 * semimajor_evolution(__orbit_power) / __semimajor;
1176  unsigned i = 1;
1177  for(; i < orbit_power_deriv.size() - 1; ++i)
1178  param_derivs[i] = semimajor_evolution(orbit_power_deriv[i]);
1179  age_deriv = semimajor_evolution(orbit_power_deriv[i]);
1180  }
1181 
1182  void BinarySystem::eccentricity_jacobian(
1183  const std::valarray<double> &orbit_power_deriv,
1184  const std::valarray<double> &orbit_angmom_gain_deriv,
1185  bool a6p5,
1186  double *param_derivs,
1187  double &age_deriv
1188  ) const
1189  {
1190  param_derivs[0] = eccentricity_evolution(__orbit_power,
1191  __orbit_angmom_gain,
1192  orbit_power_deriv[0],
1193  orbit_angmom_gain_deriv[0],
1194  true);
1195  if(a6p5) param_derivs[0] /= 6.5 * std::pow(__semimajor, 5.5);
1196  param_derivs[1] = eccentricity_evolution(__orbit_power,
1197  __orbit_angmom_gain,
1198  orbit_power_deriv[1],
1199  orbit_angmom_gain_deriv[1],
1200  false);
1201  unsigned i = 2;
1202  for(; i < orbit_power_deriv.size() - 1; ++i)
1203  param_derivs[i] = eccentricity_evolution(orbit_power_deriv[i],
1204  orbit_angmom_gain_deriv[i]);
1205  age_deriv=eccentricity_evolution(orbit_power_deriv[i],
1206  orbit_angmom_gain_deriv[i]);
1207  }
1208 
1209  void BinarySystem::angle_evolution_age_deriv(DissipatingBody &body,
1210  unsigned zone_ind,
1211  double sin_inc,
1212  double cos_inc,
1213  unsigned locked_zone_ind,
1214  double &inclination,
1215  double &periapsis) const
1216  {
1217  double dR1_dt = __body1.radius(1),
1218  dR2_dt = __body2.radius(1);
1219  DissipatingZone &zone = body.zone(zone_ind);
1220  Eigen::Vector3d orbit_torque_age_deriv = (
1221  __body1.tidal_orbit_torque(zone, Dissipation::AGE)
1222  +
1223  __body2.tidal_orbit_torque(zone, Dissipation::AGE)
1224  +
1225  __body1.tidal_orbit_torque(zone, Dissipation::RADIUS) * dR1_dt
1226  +
1227  __body2.tidal_orbit_torque(zone, Dissipation::RADIUS) * dR2_dt
1228  );
1229  DissipatingBody *other_body = &__body1;
1230  unsigned body_zone_ind = 0,
1231  num_zones = number_zones();
1232  for(
1233  unsigned other_zone_ind = 0;
1234  other_zone_ind<num_zones;
1235  ++other_zone_ind
1236  ) {
1237  orbit_torque_age_deriv += (
1238  other_body->tidal_orbit_torque(
1239  zone,
1240  Dissipation::MOMENT_OF_INERTIA,
1241  body_zone_ind,
1242  __above_lock_fractions_inertia_deriv[other_zone_ind]
1243  )
1244  *
1245  other_body->zone(body_zone_ind).moment_of_inertia(1)
1246  );
1247  ++body_zone_ind;
1248  if(body_zone_ind == other_body->number_zones()) {
1249  body_zone_ind = 0;
1250  other_body = &__body2;
1251  }
1252  }
1253  periapsis = (orbit_torque_age_deriv[1]
1254  *
1255  cos_inc
1256  /
1257  (__orbital_angmom * sin_inc));
1258  inclination = (orbit_torque_age_deriv[0] * cos_inc
1259  -
1260  orbit_torque_age_deriv[2] * sin_inc) / __orbital_angmom;
1261  Eigen::Vector3d to_add = -body.nontidal_torque(zone_ind,
1262  Dissipation::AGE);
1263  if(zone.locked()) {
1264  double above_frac =
1265  __above_lock_fractions[Dissipation::NO_DERIV][locked_zone_ind];
1266  double above_frac_deriv =
1267  __above_lock_fractions[Dissipation::AGE][locked_zone_ind];
1268  to_add -= (
1269  body.tidal_torque(zone_ind, true, Dissipation::AGE)
1270  *
1271  above_frac
1272  +
1273  body.tidal_torque(zone_ind, false, Dissipation::AGE)
1274  *
1275  (1.0 - above_frac)
1276  +
1277  above_frac_deriv*(body.tidal_torque(zone_ind, true)
1278  -
1279  body.tidal_torque(zone_ind, false))
1280  );
1281  } else to_add += body.tidal_torque(zone_ind, false, Dissipation::AGE);
1282  to_add /= zone.angular_momentum();
1283  inclination += to_add[0];
1284  periapsis -= to_add[1] / sin_inc;
1285  }
1286 
1287  void BinarySystem::angle_evolution_orbit_deriv(Dissipation::QuantityEntry entry,
1288  double angmom_deriv,
1289  DissipatingBody &body,
1290  unsigned zone_ind,
1291  double sin_inc,
1292  double cos_inc,
1293  unsigned locked_zone_ind,
1294  double &inclination,
1295  double &periapsis) const
1296  {
1297  assert(entry == Dissipation::SEMIMAJOR
1298  ||
1299  entry == Dissipation::ECCENTRICITY);
1300 
1301  DissipatingZone &zone = body.zone(zone_ind);
1302  Eigen::Vector3d orbit_torque = (__body1.tidal_orbit_torque(zone)
1303  +
1304  __body2.tidal_orbit_torque(zone)),
1305  orbit_torque_deriv = (
1306  __body1.tidal_orbit_torque(zone, entry)
1307  +
1308  __body2.tidal_orbit_torque(zone, entry)
1309  );
1310  inclination = (orbit_torque_deriv[0] * cos_inc
1311  -
1312  orbit_torque_deriv[2] * sin_inc
1313  -
1314  (orbit_torque[0] * cos_inc - orbit_torque[2] * sin_inc)
1315  *
1316  angmom_deriv
1317  /
1318  __orbital_angmom) / __orbital_angmom;
1319  periapsis = (
1320  (
1321  orbit_torque[1] * angmom_deriv / __orbital_angmom
1322  -
1323  orbit_torque_deriv[1]
1324  )
1325  *
1326  cos_inc
1327  /
1328  (__orbital_angmom * sin_inc)
1329  );
1330  Eigen::Vector3d to_add =- body.nontidal_torque(zone_ind, entry);
1331  if(zone.locked()) {
1332  double
1333  above_frac =
1334  __above_lock_fractions[Dissipation::NO_DERIV][locked_zone_ind],
1335  above_frac_deriv = __above_lock_fractions[entry][locked_zone_ind];
1336  to_add -= (
1337  body.tidal_torque(zone_ind, true, entry) * above_frac
1338  +
1339  body.tidal_torque(zone_ind, false, entry)
1340  *
1341  (1.0 - above_frac)
1342  +
1343  above_frac_deriv
1344  *
1345  (
1346  body.tidal_torque(zone_ind, true)
1347  -
1348  body.tidal_torque(zone_ind, false)
1349  )
1350  );
1351  } else to_add -= body.tidal_torque(zone_ind, false, entry);
1352  to_add /= zone.angular_momentum();
1353  inclination += to_add[0];
1354  periapsis -= to_add[1] / sin_inc;
1355  }
1356 
1357  void BinarySystem::fill_orbit_torque_deriv(
1358  Dissipation::QuantityEntry entry,
1359  DissipatingBody &body,
1360  unsigned zone_ind,
1361  std::valarray<Eigen::Vector3d> &orbit_torque_deriv
1362  ) const
1363  {
1364  assert(entry == Dissipation::INCLINATION
1365  ||
1366  entry == Dissipation::PERIAPSIS
1367  ||
1368  entry == Dissipation::MOMENT_OF_INERTIA
1369  ||
1370  entry == Dissipation::SPIN_ANGMOM);
1371  assert(orbit_torque_deriv.size()
1372  ==
1373  __body1.number_zones() + __body2.number_zones());
1374 
1375  unsigned body_deriv_zone_ind = 0,
1376  num_zones = orbit_torque_deriv.size();
1377  DissipatingBody *deriv_body = &__body1;
1378  const std::valarray<Eigen::VectorXd> *above_frac_deriv;
1379  if(entry == Dissipation::INCLINATION)
1380  above_frac_deriv = &__above_lock_fractions_inclination_deriv;
1381  else if(entry == Dissipation::PERIAPSIS)
1382  above_frac_deriv = &__above_lock_fractions_periapsis_deriv;
1383  else if(entry == Dissipation::MOMENT_OF_INERTIA)
1384  above_frac_deriv = &__above_lock_fractions_inertia_deriv;
1385  else
1386  above_frac_deriv = &__above_lock_fractions_angmom_deriv;
1387 
1388  DissipatingZone &zone = body.zone(zone_ind);
1389  for(
1390  unsigned deriv_zone_ind = 0;
1391  deriv_zone_ind < num_zones;
1392  ++deriv_zone_ind
1393  ) {
1394  orbit_torque_deriv[deriv_zone_ind] =
1395  deriv_body->tidal_orbit_torque(
1396  zone,
1397  entry,
1398  body_deriv_zone_ind,
1399  (*above_frac_deriv)[deriv_zone_ind]
1400  );
1401  if(
1402  (
1403  entry == Dissipation::INCLINATION
1404  ||
1405  entry == Dissipation::PERIAPSIS
1406  )
1407  &&
1408  deriv_body == &body
1409  &&
1410  zone_ind == body_deriv_zone_ind
1411  ) {
1412  DissipatingBody &other_body = (&body == &__body1
1413  ? __body2
1414  : __body1);
1415  orbit_torque_deriv[deriv_zone_ind] +=
1416  zone_to_zone_transform(other_body.zone(0),
1417  zone,
1418  other_body.tidal_orbit_torque(),
1419  entry,
1420  false);
1421  }
1422  ++body_deriv_zone_ind;
1423  if(body_deriv_zone_ind == deriv_body->number_zones()) {
1424  body_deriv_zone_ind = 0;
1425  deriv_body = &__body2;
1426  }
1427  }
1428  }
1429 
1430  void BinarySystem::fill_zone_torque_deriv(
1431  Dissipation::QuantityEntry entry,
1432  DissipatingBody &body,
1433  unsigned zone_ind,
1434  std::valarray<Eigen::Vector3d> &zone_torque_deriv
1435  ) const
1436  {
1437 #ifndef NDEBUG
1438  if(body.zone(zone_ind).locked())
1439  assert(zone_torque_deriv.size() == 4);
1440  else
1441  assert(zone_torque_deriv.size() == 3);
1442 #endif
1443 
1444  zone_torque_deriv[0] = (zone_ind == 0
1445  ? Eigen::Vector3d::Zero()
1446  : body.nontidal_torque(zone_ind, entry, -1));
1447  Eigen::Vector3d nontidal_torque = body.nontidal_torque(zone_ind,
1448  entry,
1449  0);
1450  zone_torque_deriv[1] = (nontidal_torque
1451  +
1452  body.tidal_torque(zone_ind, false, entry));
1453  zone_torque_deriv[2] = (zone_ind < body.number_zones() - 1
1454  ? body.nontidal_torque(zone_ind, entry, 1)
1455  : Eigen::Vector3d::Zero());
1456  if(body.zone(zone_ind).locked())
1457  zone_torque_deriv[3] = (nontidal_torque
1458  +
1459  body.tidal_torque(zone_ind, true, entry));
1460  }
1461 
1462  void BinarySystem::inclination_evolution_zone_derivs(
1463  Dissipation::QuantityEntry entry,
1464  DissipatingBody &body,
1465  unsigned zone_ind,
1466  double zone_x_torque_above,
1467  double zone_x_torque_below,
1468  const std::valarray<Eigen::Vector3d> &zone_torque_deriv,
1469  const Eigen::Vector3d &orbit_torque,
1470  const std::valarray<Eigen::Vector3d> &orbit_torque_deriv,
1471  const std::valarray<Eigen::VectorXd> &above_frac_deriv,
1472  double sin_inc,
1473  double cos_inc,
1474  unsigned locked_zone_ind,
1475  double *result
1476  ) const
1477  {
1478  assert(entry == Dissipation::INCLINATION
1479  ||
1480  entry == Dissipation::PERIAPSIS
1481  ||
1482  entry == Dissipation::SPIN_ANGMOM);
1483  assert(orbit_torque_deriv.size()
1484  ==
1485  __body1.number_zones() + __body2.number_zones());
1486 
1487  DissipatingZone &zone = body.zone(zone_ind);
1488  unsigned num_zones = orbit_torque_deriv.size();
1489  unsigned global_zone_ind = zone_ind;
1490  if(&body == &__body2) global_zone_ind += __body1.number_zones();
1491 
1492  double above_frac =
1493  __above_lock_fractions[Dissipation::NO_DERIV][locked_zone_ind];
1494 
1495  for(
1496  unsigned deriv_zone_ind = (entry == Dissipation::PERIAPSIS ? 1 : 0);
1497  deriv_zone_ind < num_zones;
1498  ++deriv_zone_ind
1499  ) {
1500  double &dest = (entry == Dissipation::PERIAPSIS
1501  ? result[deriv_zone_ind - 1]
1502  : result[deriv_zone_ind]);
1503  dest = (
1504  orbit_torque_deriv[deriv_zone_ind][0] * cos_inc
1505  -
1506  orbit_torque_deriv[deriv_zone_ind][2] * sin_inc
1507  ) / __orbital_angmom;
1508  if(zone.locked())
1509  dest -= (above_frac_deriv[deriv_zone_ind][locked_zone_ind]
1510  *
1511  (zone_x_torque_above - zone_x_torque_below)
1512  /
1513  zone.angular_momentum());
1514  if(std::abs(static_cast<int>(deriv_zone_ind - global_zone_ind)) <= 1) {
1515  double torque_deriv=
1516  zone_torque_deriv[deriv_zone_ind + 1 - global_zone_ind][0];
1517  if(zone.locked() && deriv_zone_ind == global_zone_ind) {
1518  dest-=(
1519  above_frac * zone_torque_deriv[3][0]
1520  +
1521  (1.0 - above_frac) * torque_deriv
1522  ) / zone.angular_momentum();
1523  } else dest -= torque_deriv / zone.angular_momentum();
1524  }
1525  }
1526  if(entry == Dissipation::INCLINATION)
1527  result[global_zone_ind] -= (
1528  orbit_torque[0] * sin_inc
1529  +
1530  orbit_torque[2] * cos_inc
1531  ) / __orbital_angmom;
1532  else if(entry == Dissipation::SPIN_ANGMOM)
1533  result[global_zone_ind] += (
1534  zone.locked()
1535  ? (above_frac * zone_x_torque_above
1536  +
1537  (1.0 - above_frac) * zone_x_torque_below)
1538  : zone_x_torque_below
1539  ) / std::pow(zone.angular_momentum(), 2);
1540  }
1541 
1542  void BinarySystem::periapsis_evolution_zone_derivs(
1543  Dissipation::QuantityEntry entry,
1544  DissipatingBody &body,
1545  unsigned zone_ind,
1546  double zone_y_torque_above,
1547  double zone_y_torque_below,
1548  const std::valarray<Eigen::Vector3d> &zone_torque_deriv,
1549  double orbit_y_torque,
1550  const std::valarray<Eigen::Vector3d> &orbit_torque_deriv,
1551  const std::valarray<Eigen::VectorXd> &above_frac_deriv,
1552  double sin_inc,
1553  double cos_inc,
1554  unsigned locked_zone_ind,
1555  double *result
1556  ) const
1557  {
1558  assert(entry == Dissipation::INCLINATION
1559  ||
1560  entry == Dissipation::PERIAPSIS
1561  ||
1562  entry == Dissipation::MOMENT_OF_INERTIA
1563  ||
1564  entry == Dissipation::SPIN_ANGMOM);
1565  assert(orbit_torque_deriv.size()
1566  ==
1567  __body1.number_zones() + __body2.number_zones());
1568 
1569  DissipatingZone &zone=body.zone(zone_ind);
1570  unsigned num_zones = orbit_torque_deriv.size();
1571  unsigned global_zone_ind = zone_ind;
1572  if(&body == &__body2) global_zone_ind += __body1.number_zones();
1573 
1574  double above_frac =
1575  __above_lock_fractions[Dissipation::NO_DERIV][locked_zone_ind];
1576 
1577  for(
1578  unsigned deriv_zone_ind = 0;
1579  deriv_zone_ind < num_zones;
1580  ++deriv_zone_ind
1581  ) {
1582  if(entry == Dissipation::PERIAPSIS && deriv_zone_ind == 0) continue;
1583  double &dest = (entry == Dissipation::PERIAPSIS
1584  ? result[deriv_zone_ind - 1]
1585  : result[deriv_zone_ind]);
1586  dest = (-orbit_torque_deriv[deriv_zone_ind][1]
1587  *
1588  cos_inc
1589  /
1590  (__orbital_angmom * sin_inc));
1591  if(zone.locked())
1592  dest += (above_frac_deriv[deriv_zone_ind][locked_zone_ind]
1593  *
1594  (zone_y_torque_above - zone_y_torque_below)
1595  /
1596  (sin_inc * zone.angular_momentum()));
1597  if(
1598  std::abs(static_cast<int>(deriv_zone_ind - global_zone_ind))
1599  <=
1600  1
1601  ) {
1602  double torque_deriv =
1603  zone_torque_deriv[deriv_zone_ind + 1 - global_zone_ind][1];
1604  if(zone.locked() && deriv_zone_ind == global_zone_ind) {
1605  dest+=(
1606  above_frac * zone_torque_deriv[3][1]
1607  +
1608  (1.0 - above_frac) * torque_deriv
1609  ) / (sin_inc * zone.angular_momentum());
1610  } else dest += (torque_deriv
1611  /
1612  (sin_inc * zone.angular_momentum()));
1613  }
1614  }
1615  double zone_torque;
1616  if(zone.locked()) zone_torque = (above_frac*zone_y_torque_above
1617  +
1618  (1.0 - above_frac)
1619  *
1620  zone_y_torque_below);
1621  else zone_torque = zone_y_torque_below;
1622  if(entry == Dissipation::INCLINATION)
1623  result[global_zone_ind] += (
1624  orbit_y_torque / (__orbital_angmom * std::pow(sin_inc, 2))
1625  -
1626  zone_torque*cos_inc
1627  /
1628  (zone.angular_momentum() * std::pow(sin_inc, 2))
1629  );
1630  else if(entry == Dissipation::SPIN_ANGMOM)
1631  result[global_zone_ind] -=
1632  zone_torque / (std::pow(zone.angular_momentum(), 2) * sin_inc);
1633  }
1634 
1635  void BinarySystem::spin_angmom_evolution_zone_derivs(
1636  Dissipation::QuantityEntry entry,
1637  DissipatingBody &body,
1638  unsigned zone_ind,
1639  double zone_z_torque_above,
1640  double zone_z_torque_below,
1641  const std::valarray<Eigen::Vector3d> &zone_torque_deriv,
1642  const std::valarray<Eigen::VectorXd> &above_frac_deriv,
1643  unsigned locked_zone_ind,
1644  double *result
1645  ) const
1646  {
1647  unsigned global_zone_ind = zone_ind,
1648  num_zones = number_zones();
1649  if(&body == &__body2) global_zone_ind += __body1.number_zones();
1650  double above_frac =
1651  __above_lock_fractions[Dissipation::NO_DERIV][locked_zone_ind];
1652  bool zone_is_locked = body.zone(zone_ind).locked();
1653  for(
1654  unsigned deriv_zone_ind = 0;
1655  deriv_zone_ind<num_zones;
1656  ++deriv_zone_ind
1657  ) {
1658  double &dest = (entry == Dissipation::PERIAPSIS
1659  ? result[deriv_zone_ind - 1]
1660  : result[deriv_zone_ind]);
1661  if(zone_is_locked)
1662  dest = (above_frac_deriv[deriv_zone_ind][locked_zone_ind]
1663  *
1664  (zone_z_torque_above - zone_z_torque_below));
1665  else dest = 0;
1666  if(std::abs(static_cast<int>(deriv_zone_ind
1667  -
1668  global_zone_ind)) <= 1) {
1669  double torque_deriv =
1670  zone_torque_deriv[deriv_zone_ind + 1 - global_zone_ind][2];
1671  if(zone_is_locked) dest += (above_frac * zone_torque_deriv[3][2]
1672  +
1673  (1.0 - above_frac) * torque_deriv);
1674  else dest += torque_deriv;
1675  }
1676  }
1677  }
1678 
1679  void BinarySystem::binary_jacobian(double *param_derivs,
1680  double *age_derivs) const
1681  {
1682  unsigned body1_zones = __body1.number_zones(),
1683  num_zones = body1_zones + __body2.number_zones(),
1684  num_locked_zones = (__body1.number_locked_zones()
1685  +
1686  __body2.number_locked_zones()),
1687  num_param = 1 + 6 * num_zones - num_locked_zones;
1688  double dangmom_da = __orbital_angmom / (2.0 * __semimajor),
1689  dangmom_de = (-__eccentricity
1690  /
1691  (1.0 - std::pow(__eccentricity, 2))
1692  *
1693  __orbital_angmom);
1694  std::valarray<double> orbit_power_deriv(num_param + 1),
1695  orbit_angmom_gain_deriv(num_param + 1);
1696  fill_orbit_power_deriv(orbit_power_deriv);
1697  fill_orbit_angmom_gain_deriv(orbit_angmom_gain_deriv);
1698  semimajor_jacobian(orbit_power_deriv,
1699  num_locked_zones == 0,
1700  param_derivs,
1701  age_derivs[0]);
1702  eccentricity_jacobian(orbit_power_deriv,
1703  orbit_angmom_gain_deriv,
1704  num_locked_zones == 0,
1705  param_derivs + num_param,
1706  age_derivs[1]);
1707  unsigned locked_zone_ind = 0;
1708  DissipatingBody *body = &__body1;
1709  unsigned body_zone_ind = 0;
1710  static Dissipation::QuantityEntry zone_deriv_list[]={
1711  Dissipation::INCLINATION,
1712  Dissipation::PERIAPSIS,
1713  Dissipation::SPIN_ANGMOM
1714  };
1715  std::valarray<double> reference_periapsis_param_deriv(num_param);
1716  double reference_periapsis_age_deriv;
1717  const std::valarray<Eigen::VectorXd> *above_frac_deriv[] = {
1718  &__above_lock_fractions_inclination_deriv,
1719  &__above_lock_fractions_periapsis_deriv,
1720  &__above_lock_fractions_angmom_deriv
1721  };
1722  for(unsigned zone_ind = 0; zone_ind < num_zones; ++zone_ind) {
1723  DissipatingZone &zone = body->zone(zone_ind);
1724  double sin_inc = std::sin(zone.inclination()),
1725  cos_inc = std::cos(zone.inclination());
1726  double &periapsis_age_deriv = (
1727  zone_ind==0
1728  ? reference_periapsis_age_deriv
1729  : age_derivs[1 + num_zones + zone_ind]
1730  );
1731  angle_evolution_age_deriv(*body,
1732  zone_ind,
1733  sin_inc,
1734  cos_inc,
1735  locked_zone_ind,
1736  age_derivs[2 + zone_ind],
1737  periapsis_age_deriv);
1738 
1739  double *periapsis_row = (
1740  zone_ind==0
1741  ? &reference_periapsis_param_deriv[0]
1742  : param_derivs+(1 + zone_ind+num_zones) * num_param
1743  );
1744 
1745  angle_evolution_orbit_deriv(Dissipation::SEMIMAJOR,
1746  dangmom_da,
1747  *body,
1748  body_zone_ind,
1749  sin_inc,
1750  cos_inc,
1751  locked_zone_ind,
1752  param_derivs[(2 + zone_ind) * num_param],
1753  periapsis_row[0]);
1754  angle_evolution_orbit_deriv(
1755  Dissipation::ECCENTRICITY,
1756  dangmom_de,
1757  *body,
1758  body_zone_ind,
1759  sin_inc, cos_inc,
1760  locked_zone_ind,
1761  param_derivs[(2 + zone_ind) * num_param + 1],
1762  periapsis_row[1]
1763  );
1764 
1765  Eigen::Vector3d zone_torque_above = body->nontidal_torque(zone_ind),
1766  zone_torque_below = zone_torque_above,
1767  orbit_torque = body->tidal_orbit_torque(zone);
1768  zone_torque_above += body->tidal_torque(zone_ind, true);
1769  zone_torque_below += body->tidal_torque(zone_ind, false);
1770  std::valarray<Eigen::Vector3d>
1771  zone_torque_deriv(zone.locked() ? 4 : 3),
1772  orbit_torque_deriv(num_zones);
1773  unsigned param_offset = (2 + zone_ind) * num_param + 2;
1774  for(unsigned deriv_ind = 0; deriv_ind < 3; ++deriv_ind) {
1775  Dissipation::QuantityEntry deriv = zone_deriv_list[deriv_ind];
1776  fill_zone_torque_deriv(deriv,
1777  *body,
1778  zone_ind,
1779  zone_torque_deriv);
1780  fill_orbit_torque_deriv(deriv,
1781  *body,zone_ind,
1782  orbit_torque_deriv);
1783  inclination_evolution_zone_derivs(deriv,
1784  *body, zone_ind,
1785  zone_torque_above[0],
1786  zone_torque_below[0],
1787  zone_torque_deriv,
1788  orbit_torque,
1789  orbit_torque_deriv,
1790  *(above_frac_deriv[deriv_ind]),
1791  sin_inc,
1792  cos_inc,
1793  locked_zone_ind,
1794  param_derivs + param_offset);
1795  periapsis_evolution_zone_derivs(deriv,
1796  *body,
1797  zone_ind,
1798  zone_torque_above[1],
1799  zone_torque_below[1],
1800  zone_torque_deriv,
1801  orbit_torque[1],
1802  orbit_torque_deriv,
1803  *(above_frac_deriv[deriv_ind]),
1804  sin_inc,
1805  cos_inc,
1806  locked_zone_ind,
1807  periapsis_row + 2);
1808  spin_angmom_evolution_zone_derivs(
1809  deriv,
1810  *body,
1811  zone_ind,
1812  zone_torque_above[2],
1813  zone_torque_below[2],
1814  zone_torque_deriv,
1815  *(above_frac_deriv[deriv_ind]),
1816  locked_zone_ind,
1817  param_derivs + param_offset + 2 * num_zones - 1
1818  );
1819  param_offset += (deriv == Dissipation::PERIAPSIS
1820  ? num_zones - 1
1821  : num_zones);
1822  }
1823 
1824  if(zone.locked()) ++locked_zone_ind;
1825  ++body_zone_ind;
1826  if(body_zone_ind == body->number_zones()) {
1827  body_zone_ind = 0;
1828  body = &__body2;
1829  }
1830  if(zone_ind != 0) {
1831  periapsis_age_deriv -= reference_periapsis_age_deriv;
1832  for(unsigned i = 0; i < num_param; ++i)
1833  periapsis_row[i] -= reference_periapsis_param_deriv[i];
1834  }
1835  }
1836 
1837  }
1838 
1839  void BinarySystem::fill_locked_surface_orbit(
1840  std::valarray<double> &orbit
1841  ) const
1842  {
1843  assert(__evolution_mode == Core::LOCKED_SURFACE_SPIN);
1844  assert(std::isnan(__semimajor));
1845  assert(std::isnan(__eccentricity));
1846 
1847  orbit.resize(__body1.number_zones() - 1);
1848  for(unsigned i = 1; i < __body1.number_zones(); ++i) {
1849  orbit[i - 1] = __body1.zone(i).angular_momentum();
1850  assert(__body1.zone(i).inclination() == 0);
1851  assert(__body1.zone(i).periapsis() == 0);
1852  }
1853  }
1854 
1855  void BinarySystem::fill_binary_orbit(std::valarray<double> &orbit) const
1856  {
1857  assert(__evolution_mode == Core::BINARY);
1858  assert(!std::isnan(__semimajor));
1859  assert(!std::isnan(__eccentricity));
1860  assert(__body1.zone(0).periapsis() == 0);
1861 
1862  orbit.resize(1 + 3 * number_zones() - number_locked_zones());
1863  if(number_locked_zones() == 0) orbit[0] = std::pow(__semimajor, 6.5);
1864  else orbit[0] = __semimajor;
1865  orbit[1] = __eccentricity;
1866  unsigned inclination_ind = 2,
1867  periapsis_ind = 2 + number_zones(),
1868  angmom_ind = periapsis_ind + number_zones() - 1;
1869  for(short body_ind = 0; body_ind < 2; ++body_ind) {
1870  DissipatingBody &body = (body_ind == 0 ? __body1 : __body2);
1871  for(
1872  unsigned zone_ind = 0;
1873  zone_ind < body.number_zones();
1874  ++zone_ind
1875  ) {
1876  DissipatingZone &zone = body.zone(zone_ind);
1877  orbit[inclination_ind++] = zone.inclination();
1878  if(body_ind || zone_ind)
1879  orbit[periapsis_ind++] = zone.periapsis();
1880  if(!zone.locked())
1881  orbit[angmom_ind++] = zone.angular_momentum();
1882  }
1883  }
1884  }
1885 
1886  void BinarySystem::fill_single_orbit(std::valarray<double> &orbit) const
1887  {
1888  assert(__evolution_mode == Core::SINGLE);
1889  assert(std::isnan(__semimajor));
1890  assert(std::isnan(__eccentricity));
1891  assert(__body1.zone(0).periapsis() == 0);
1892  assert(__body1.zone(0).inclination() == 0);
1893  assert(__body1.number_locked_zones() == 0);
1894 
1895  orbit.resize(3 * __body1.number_zones() - 2);
1896  unsigned inclination_ind = 0,
1897  periapsis_ind = __body1.number_zones() - 1,
1898  angmom_ind = periapsis_ind + __body1.number_zones() - 1;
1899  for(
1900  unsigned zone_ind = 0;
1901  zone_ind < __body1.number_zones();
1902  ++zone_ind
1903  ) {
1904  DissipatingZone &zone = __body1.zone(zone_ind);
1905  if(zone_ind) orbit[inclination_ind++] = zone.inclination();
1906  if(zone_ind) orbit[periapsis_ind++] = zone.periapsis();
1907  assert(!zone.locked());
1908  orbit[angmom_ind++] = zone.angular_momentum();
1909  }
1910  }
1911 
1912  int BinarySystem::configure(bool initialize,
1913  double age,
1914  double semimajor,
1915  double eccentricity,
1916  const double *spin_angmom,
1917  const double *inclination,
1918  const double *periapsis,
1919  Core::EvolModeType evolution_mode)
1920  {
1921 #ifndef NDEBUG
1922  if(initialize)
1923  std::cerr << "Initializing BinarySystem." << std::endl;
1924  if(evolution_mode != Core::BINARY) {
1925  assert(std::isnan(semimajor));
1926  assert(std::isnan(eccentricity));
1927  }
1928  if(evolution_mode == Core::LOCKED_SURFACE_SPIN) {
1929  assert(inclination == NULL);
1930  assert(periapsis == NULL);
1931  }
1932 
1933  std::cerr << "Configuring binary with a = " << semimajor
1934  << ", e = " << eccentricity
1935  << " in " << evolution_mode << " mode"
1936  << std::endl;
1937 #endif
1938 
1939  if(
1940  evolution_mode == Core::BINARY
1941  &&
1942  !(
1943  semimajor > 0.0
1944  &&
1945  eccentricity >= 0.0
1946  &&
1947  eccentricity < 1.0
1948  )
1949  )
1950  return GSL_EDOM;
1951 
1952  __evolution_mode = evolution_mode;
1953  double m1 = __body1.mass(),
1954  m2 = __body2.mass();
1955  __age = age;
1956  __semimajor = semimajor;
1957  __eccentricity = eccentricity;
1958  __orbital_energy = Core::orbital_energy(m1, m2, semimajor);
1959  __orbital_angmom = Core::orbital_angular_momentum(m1,
1960  m2,
1961  semimajor,
1962  eccentricity);
1963 #ifndef NDEBUG
1964  std::cerr << "Configuring primary." << std::endl;
1965 #endif
1966  __body1.configure(initialize,
1967  age,
1968  m2,
1969  semimajor,
1970  eccentricity,
1971  spin_angmom,
1972  inclination,
1973  periapsis,
1974  evolution_mode == Core::LOCKED_SURFACE_SPIN,
1975  evolution_mode != Core::BINARY,
1976  true);
1977 
1978 
1979  if(evolution_mode == Core::BINARY) {
1980  unsigned offset = __body1.number_zones();
1981 #ifndef NDEBUG
1982  std::cerr << "Configuring secondary." << std::endl;
1983 #endif
1984  __body2.configure(initialize,
1985  age,
1986  m1,
1987  semimajor,
1988  eccentricity,
1989  spin_angmom + offset - __body1.number_locked_zones(),
1990  inclination + offset,
1991  periapsis + offset - 1);
1992  find_locked_zones();
1993  update_above_lock_fractions();
1994  } else
1995  find_locked_zones();
1996 
1997  fill_orbit_torque_and_power();
1998 #ifndef NDEBUG
1999  // if(evolution_mode == BINARY)
2000  // assert(__semimajor > minimum_semimajor());
2001 #endif
2002  return 0;
2003  }
2004 
2005  int BinarySystem::configure(bool initialize,
2006  double age,
2007  const double *parameters,
2008  Core::EvolModeType evolution_mode)
2009  {
2010  __evolution_mode = evolution_mode;
2011  double semimajor, eccentricity;
2012  const double *spin_angmom, *inclination, *periapsis;
2013  unsigned num_zones = number_zones();
2014  if(evolution_mode == Core::BINARY) {
2015  if(parameters[0] < 0) {
2016 #ifndef NDEBUG
2017  std::cerr << "At t = " << age << " param: ";
2018  for(
2019  unsigned i = 0;
2020  i < 3 * num_zones + 1 - number_locked_zones();
2021  ++i
2022  ) {
2023  if(i) std::cerr << ", ";
2024  std::cerr << parameters[i];
2025  }
2026  std::cerr << std::endl;
2027 #endif
2028  return GSL_EDOM;
2029  }
2030  if(__body1.number_locked_zones() || __body2.number_locked_zones()) {
2031  semimajor = parameters[0];
2032  } else
2033  semimajor = std::pow(parameters[0], 1.0 / 6.5);
2034  eccentricity = parameters[1];
2035  inclination = parameters+2;
2036  } else {
2037  semimajor = eccentricity=Core::NaN;
2038  if(evolution_mode == Core::SINGLE) inclination = parameters;
2039  else inclination = NULL;
2040  }
2041 
2042  if(evolution_mode == Core::LOCKED_SURFACE_SPIN) {
2043  periapsis = NULL;
2044  spin_angmom = parameters;
2045  } else {
2046  assert(inclination != NULL);
2047 
2048  periapsis = inclination+num_zones;
2049  if(evolution_mode == Core::SINGLE) --periapsis;
2050  spin_angmom = periapsis + num_zones - 1;
2051  }
2052  return configure(initialize,
2053  age,
2054  semimajor,
2055  eccentricity,
2056  spin_angmom,
2057  inclination,
2058  periapsis,
2059  evolution_mode);
2060  }
2061 
2062  Core::EvolModeType BinarySystem::fill_orbit(
2063  std::valarray<double> &orbit
2064  ) const
2065  {
2066  if(__evolution_mode == Core::LOCKED_SURFACE_SPIN)
2067  fill_locked_surface_orbit(orbit);
2068  else if(__evolution_mode == Core::BINARY) fill_binary_orbit(orbit);
2069  else fill_single_orbit(orbit);
2070  return __evolution_mode;
2071  }
2072 
2073  double BinarySystem::above_lock_fraction(unsigned locked_zone_index,
2074  Dissipation::QuantityEntry entry,
2075  unsigned deriv_zone_index,
2076  bool secondary_radius)
2077  {
2078  if(zone_specific(entry)) {
2079  assert(entry == Dissipation::INCLINATION
2080  ||
2081  entry == Dissipation::PERIAPSIS
2082  ||
2083  entry == Dissipation::SPIN_ANGMOM
2084  ||
2085  entry == Dissipation::MOMENT_OF_INERTIA);
2086 
2087  if(entry == Dissipation::INCLINATION)
2088  return __above_lock_fractions_inclination_deriv[deriv_zone_index]
2089  [locked_zone_index];
2090  else if(entry == Dissipation::PERIAPSIS)
2091  return __above_lock_fractions_periapsis_deriv[deriv_zone_index]
2092  [locked_zone_index];
2093  else if(entry == Dissipation::SPIN_ANGMOM)
2094  return __above_lock_fractions_angmom_deriv[deriv_zone_index]
2095  [locked_zone_index];
2096  else return __above_lock_fractions_inertia_deriv[deriv_zone_index]
2097  [locked_zone_index];
2098  } else {
2099  if(entry == Dissipation::RADIUS && secondary_radius)
2100  return
2101  __above_lock_fractions_body2_radius_deriv[locked_zone_index];
2102  else return __above_lock_fractions[entry][locked_zone_index];
2103  }
2104  }
2105 
2106  int BinarySystem::differential_equations(double age,
2107  const double *parameters,
2108  Core::EvolModeType evolution_mode,
2109  double *differential_equations,
2110  bool expansion_error)
2111  {
2112 #ifndef NDEBUG
2113  if(!expansion_error)
2114  std::cerr << "Finding differential equations at t = " << age
2115  << " in " << evolution_mode
2116  << " mode, with orbit[0] = " << parameters[0]
2117  << std::endl;
2118 #endif
2119  int status = configure(false, age, parameters, evolution_mode);
2120  if(status != GSL_SUCCESS) return status;
2121  if(!expansion_error)
2122  __semimajor_rate = __eccentricity_rate = Core::NaN;
2123  switch(evolution_mode) {
2124  case Core::LOCKED_SURFACE_SPIN :
2125  return locked_surface_differential_equations(
2126  differential_equations,
2127  expansion_error
2128  );
2129  case Core::SINGLE :
2130  return single_body_differential_equations(
2131  differential_equations,
2132  expansion_error
2133  );
2134  case Core::BINARY :
2135  return binary_differential_equations(differential_equations,
2136  expansion_error);
2137  default :
2138  throw Core::Error::BadFunctionArguments(
2139  "Evolution mode other than LOCKED_SURFACE_SPIN, SINGLE or "
2140  "BINARY encountered in "
2141  "BinarySystem::differential_equations!"
2142  );
2143  }
2144  }
2145 
2146  int BinarySystem::jacobian(double age,
2147  const double *parameters,
2148  Core::EvolModeType evolution_mode,
2149  double *param_derivs,
2150  double *age_derivs)
2151  {
2152  configure(false, age, parameters, evolution_mode);
2153  switch(evolution_mode) {
2154  case Core::LOCKED_SURFACE_SPIN :
2155  locked_surface_jacobian(param_derivs, age_derivs);
2156  return 0;
2157  case Core::SINGLE : single_body_jacobian(param_derivs, age_derivs);
2158  return 0;
2159  case Core::BINARY : binary_jacobian(param_derivs, age_derivs);
2160  return 0;
2161  default : throw Core::Error::BadFunctionArguments(
2162  "Evolution mode other than LOCKED_SURFACE_SPIN, "
2163  "SINGLE or BINARY encountered in "
2164  "BinarySystem::jacobian!"
2165  );
2166  }
2167  }
2168 
2169  void BinarySystem::check_for_lock(int orbital_freq_mult,
2170  int spin_freq_mult,
2171  unsigned short body_index,
2172  unsigned zone_index,
2173  short direction)
2174  {
2175  DissipatingBody &body = (body_index ? __body2 : __body1);
2176 
2177  assert(body_index <= 1);
2178  assert(zone_index < body.number_zones());
2179  assert(spin_freq_mult);
2180 
2181  double original_angmom = body.zone(zone_index).angular_momentum();
2182  body.lock_zone_spin(zone_index, orbital_freq_mult, spin_freq_mult);
2183  unsigned num_zones = number_zones(),
2184  num_locked_zones = number_locked_zones(),
2185  locked_zone_ind = 0;
2186  std::vector<double> spin_angmom(num_zones - num_locked_zones),
2187  inclinations(num_zones),
2188  periapses(num_zones);
2189  for(unsigned zone_ind = 0; zone_ind < num_zones; ++zone_ind) {
2190  DissipatingZone &zone=(zone_ind < __body1.number_zones()
2191  ? __body1.zone(zone_ind)
2192  : __body2.zone(zone_ind
2193  -
2194  __body1.number_zones()));
2195  if(zone.locked()) ++locked_zone_ind;
2196  else spin_angmom[zone_ind-locked_zone_ind] = (
2197  zone.angular_momentum()
2198  );
2199  inclinations[zone_ind] = zone.inclination();
2200  if(zone_ind) periapses[zone_ind - 1] = zone.periapsis();
2201  }
2202  assert(locked_zone_ind == num_locked_zones);
2203  configure(false,
2204  __age,
2205  __semimajor,
2206  __eccentricity,
2207  &(spin_angmom[0]),
2208  &(inclinations[0]),
2209  &(periapses[0]),
2210  Core::BINARY);
2211  DissipatingZone &locked_zone = body.zone(zone_index);
2212  double above_lock_fraction = __above_lock_fractions
2213  [Dissipation::NO_DERIV]
2214  [locked_zone.locked_zone_index()];
2215 #ifndef NDEBUG
2216  std::cerr << "Holding lock requires above lock fraction of: "
2217  << above_lock_fraction
2218  << std::endl;
2219 #endif
2220  if(above_lock_fraction > 0 && above_lock_fraction < 1) return;
2221  std::vector<double>::iterator check_zone_dest = (
2222  spin_angmom.begin()
2223  +
2224  (
2225  zone_index
2226  +
2227  (body_index ? __body1.number_zones() : 0)
2228  -
2229  locked_zone.locked_zone_index()
2230  )
2231  );
2232  spin_angmom.insert(check_zone_dest, original_angmom);
2233  if(direction == 0)
2234  throw Core::Error::Runtime(
2235  "Crossing spin-orbit synchronization with unspecified direction!"
2236  );
2237  body.unlock_zone_spin(zone_index, direction);
2238  configure(false,
2239  __age,
2240  __semimajor,
2241  __eccentricity,
2242  &(spin_angmom[0]),
2243  &(inclinations[0]),
2244  &(periapses[0]),
2245  Core::BINARY);
2246  }
2247 
2248  double BinarySystem::minimum_separation(bool deriv) const
2249  {
2250  double rroche = (
2251  __body2.radius()
2252  ?(
2253  2.44
2254  *
2255  __body2.radius()
2256  *
2257  std::pow(__body1.mass() / __body2.mass(), 1.0 / 3.0)
2258  )
2259  :(
2260  0.0
2261  )
2262  );
2263  if(rroche > __body1.radius()) {
2264  if(deriv) return rroche * __body2.radius(1) / __body2.radius();
2265  else return rroche;
2266  } else {
2267  return __body1.radius(deriv ? 1 : 0);
2268  }
2269  }
2270 
2271  void BinarySystem::secondary_died()
2272  {
2273 #ifndef NDEBUG
2274  std::cerr << "Handling secondary death!" << std::endl;
2275 #endif
2276  unsigned num_zones = __body1.number_zones();
2277  std::valarray<double> spin_angmom(num_zones),
2278  inclination(num_zones - 1),
2279  periapsis(num_zones - 1);
2280  DissipatingZone &old_surface_zone = __body1.zone(0);
2281  double old_surface_inclination = old_surface_zone.inclination(),
2282  sin_inc = std::sin(old_surface_inclination),
2283  cos_inc = std::cos(old_surface_inclination),
2284  old_surface_angmom = old_surface_zone.angular_momentum(),
2285  angmom = std::sqrt(
2286  std::pow(old_surface_angmom + __orbital_angmom * cos_inc, 2)
2287  +
2288  std::pow(__orbital_angmom * sin_inc, 2)
2289  ),
2290  new_surface_inclination=std::atan2(__orbital_angmom*sin_inc,
2291  old_surface_angmom
2292  +
2293  __orbital_angmom*cos_inc);
2294 
2295  spin_angmom[0] = angmom;
2296  // assert(num_zones == 2);
2297  if(old_surface_zone.locked()) __body1.unlock_zone_spin(0, 1);
2298  for(unsigned zone_ind = 1; zone_ind < num_zones; ++zone_ind) {
2299  DissipatingZone &zone = __body1.zone(zone_ind);
2300  assert(zone.periapsis() == 0);
2301  inclination[zone_ind - 1]=std::abs(zone.inclination()
2302  -
2303  old_surface_inclination
2304  +
2305  new_surface_inclination);
2306  periapsis[zone_ind - 1] = 0;
2307  spin_angmom[zone_ind] = zone.angular_momentum();
2308  if(zone.locked()) __body1.unlock_zone_spin(zone_ind, 1);
2309  }
2310  configure(false,
2311  __age,
2312  Core::NaN,
2313  Core::NaN,
2314  &spin_angmom[0],
2315  &inclination[0],
2316  &periapsis[0],
2317  Core::SINGLE);
2318  __body1.spin_jumped();
2319  }
2320 
2321  void BinarySystem::release_lock(unsigned locked_zone_index, short direction)
2322  {
2323  DissipatingBody *body;
2324  if(locked_zone_index < __body1.number_locked_zones()) {
2325  body = &__body1;
2326  for(
2327  unsigned zone_ind = 0;
2328  zone_ind < __body2.number_locked_zones();
2329  ++zone_ind
2330  )
2331  if(__body2.zone(zone_ind).locked())
2332  --__body2.zone(zone_ind).locked_zone_index();
2333 
2334  } else {
2335  body = &__body2;
2336  locked_zone_index -= __body1.number_locked_zones();
2337  }
2338  unsigned zone_ind = 0;
2339  while(true){
2340  if(body->zone(zone_ind).locked()) {
2341  if(locked_zone_index == 0) break;
2342  else --locked_zone_index;
2343  }
2344  assert(zone_ind < body->number_zones());
2345  ++zone_ind;
2346  }
2347  body->unlock_zone_spin(zone_ind, direction);
2348  for(; zone_ind < body->number_zones(); ++zone_ind)
2349  if(body->zone(zone_ind).locked())
2350  --body->zone(zone_ind).locked_zone_index();
2351  }
2352 
2353  void BinarySystem::add_to_evolution()
2354  {
2355  __semimajor_evolution.push_back(__semimajor);
2356  __eccentricity_evolution.push_back(__eccentricity);
2357  __semimajor_rate_evolution.push_back(__semimajor_rate);
2358  __eccentricity_rate_evolution.push_back(__eccentricity_rate);
2359  __body1.add_to_evolution();
2360  __body2.add_to_evolution();
2361  }
2362 
2363  void BinarySystem::reset_evolution()
2364  {
2365  __semimajor_evolution.clear();
2366  __eccentricity_evolution.clear();
2367  __semimajor_rate_evolution.clear();
2368  __eccentricity_rate_evolution.clear();
2369  __body1.reset_evolution();
2370  __body2.reset_evolution();
2371  }
2372 
2373  void BinarySystem::rewind_evolution(unsigned nsteps)
2374  {
2375  for(unsigned i = 0; i < nsteps; ++i) {
2376  __semimajor_evolution.pop_back();
2377  __eccentricity_evolution.pop_back();
2378  __semimajor_rate_evolution.pop_back();
2379  __eccentricity_rate_evolution.pop_back();
2380  }
2381  __body1.rewind_evolution(nsteps);
2382  __body2.rewind_evolution(nsteps);
2383  }
2384 
2385  CombinedStoppingCondition *BinarySystem::stopping_conditions()
2386  {
2387  CombinedStoppingCondition *result = new CombinedStoppingCondition();
2388  if(__evolution_mode == Core::BINARY)
2389  (*result) |= new SecondaryDeathCondition(*this);
2390  (*result) |= __body1.stopping_conditions(*this, true);
2391  if(__evolution_mode == Core::BINARY)
2392  (*result) |= __body2.stopping_conditions(*this, false);
2393  return result;
2394  }
2395 
2396  void BinarySystem::reached_critical_age(double age)
2397  {
2398  __body1.reached_critical_age(age);
2399  __body2.reached_critical_age(age);
2400  }
2401 
2402  double BinarySystem::next_stop_age() const
2403  {
2404  return std::min(__body1.next_stop_age(),
2405  __body2.next_stop_age());
2406  }
2407 
2408  unsigned BinarySystem::eccentricity_order() const
2409  {
2410 #ifndef NDEBUG
2411  int result = -1;
2412 #endif
2413  DissipatingBody *body = &__body1;
2414  while(true) {
2415  for(
2416  unsigned zone_ind = 0;
2417  zone_ind < body->number_zones();
2418  ++zone_ind
2419  ) {
2420  DissipatingZone &zone = body->zone(zone_ind);
2421  if(zone.dissipative()) {
2422 #ifndef NDEBUG
2423  if(result < 1)
2424  result = zone.eccentricity_order();
2425  else
2426  assert(static_cast<unsigned>(result)
2427  ==
2428  zone.eccentricity_order());
2429 #else
2430  return zone.eccentricity_order();
2431 #endif
2432  }
2433  }
2434  if(body == &__body2) break;
2435  body = &__body2;
2436  }
2437 #ifndef NDEBUG
2438  return result;
2439 #endif
2440  return 0;
2441  }
2442 
2443 } //End Evolve namespace.
Evolve::BinarySystem::__orbit_angmom_gain_expansion_error
double __orbit_angmom_gain_expansion_error
Estimate of the error in __orbit_angmom_gain due to truncating the tidal potential eccentricity expan...
Definition: BinarySystem.h:80
Evolve::BinarySystem::add_to_evolution
virtual void add_to_evolution()
Appends the state defined by last configure(), to the evolution.
Definition: BinarySystem.cpp:2353
Evolve::BinarySystem::__age
double __age
The present age of the stellar system in Gyrs.
Definition: BinarySystem.h:80
Evolve::BinarySystem::fill_single_body_jacobian
void fill_single_body_jacobian(double *inclination_param_derivs, double *periapsis_param_derivs, double *angmom_param_derivs, double *inclination_age_derivs, double *periapsis_age_derivs, double *angmom_age_derivs) const
Fills the jacobian for a system consisting of one isolated body.
Definition: BinarySystem.cpp:229
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::BinarySystem::differential_equations
int differential_equations(double age, const double *parameters, Core::EvolModeType evolution_mode, double *differential_equations, bool expansion_error=false)
The differential equation and jacobian for the evolution of the system.
Definition: BinarySystem.cpp:2106
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
Evolve::BinarySystem::__above_lock_fractions_body2_radius_deriv
Eigen::VectorXd __above_lock_fractions_body2_radius_deriv
The derivative of the above lock fractions w.r.t. to the radius of the secondardy.
Definition: BinarySystem.h:159
Core::Error::BadFunctionArguments
Function arguments do not satisfy some requirement.
Definition: Error.h:73
RADIUS
RADIUS
The derivative w.r.t. the radius of the body in .
Definition: DissipationQuantities.h:99
Evolve::BinarySystem::eccentricity_evolution
const std::list< double > & eccentricity_evolution() const
The tabulated evolution of the eccentricity so far.
Definition: BinarySystem.h:1052
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::BinarySystem::__orbit_torque
Eigen::Vector3d __orbit_torque
The torque on the orbit in the coordinate system of the outermost zone of the first body...
Definition: BinarySystem.h:119
Evolve::BinarySystem::find_locked_zones
void find_locked_zones()
Fills the __locked_zones list.
Definition: BinarySystem.cpp:13
Evolve::DissipatingZone::angular_momentum
double angular_momentum() const
The angular momentum of the given zone in .
Definition: DissipatingZone.h:587
Evolve::BinarySystem::__semimajor_rate
double __semimajor_rate
The current rate at which the semiamjor axis is changing.
Definition: BinarySystem.h:111
Evolve::DissipatingBody::number_locked_zones
unsigned number_locked_zones() const
The number of zones currently in a spin-orbit lock.
Definition: DissipatingBody.h:339
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::BinarySystem::single_body_differential_equations
int single_body_differential_equations(double *evolution_rates, bool expansion_error) const
Differential equations for the rotation of the zones of body 0 if no other body is present...
Definition: BinarySystem.cpp:163
Evolve::BinarySystem::__eccentricity_evolution
std::list< double > __eccentricity_evolution
The evolution of the eccentricity recorded by add_to_evolution() so far.
Definition: BinarySystem.h:64
SPIN_FREQUENCY
SPIN_FREQUENCY
The derivative w.r.t. the spin frequency of a dissipating zone.
Definition: DissipationQuantities.h:63
Evolve::BinarySystem::__semimajor_evolution
std::list< double > __semimajor_evolution
The evolution of the semimajor axis recorded by add_to_evolution() so far.
Definition: BinarySystem.h:64
Evolve::DissipatingBody
A base class for any body contributing to tidal dissipation.
Definition: DissipatingBody.h:39
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::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
Evolve::BinarySystem::fill_above_lock_fractions_deriv
void fill_above_lock_fractions_deriv()
Fills the __above_lock_fractinos_*_deriv members.
Definition: BinarySystem.cpp:779
Evolve::ZoneOrientation::periapsis
double periapsis(bool evolution_rate=false) const
The argument of periapsis of this zone minus the reference zone&#39;s.
Definition: ZoneOrientation.h:59
Evolve::BinarySystem::binary_differential_equations
int binary_differential_equations(double *differential_equations, bool expansion_error) const
The differential equations for a system with both bodies present.
Definition: BinarySystem.cpp:895
Evolve::DissipatingBody::zone
virtual const DissipatingZone & zone(unsigned zone_index) const =0
A modifiable reference to one of the body&#39;s zones.
Evolve::BinarySystem::__semimajor
double __semimajor
The current semimajor axis.
Definition: BinarySystem.h:80
Evolve::BinarySystem::check_for_lock
void check_for_lock(int orbital_freq_mult, int spin_freq_mult, unsigned short body_index, unsigned zone_index, short direction)
Check if a spin-orbit lock can be held and updates the system as necessary to calculate subsequent ev...
Definition: BinarySystem.cpp:2169
Evolve::BinarySystem::__above_lock_fractions
std::valarray< Eigen::VectorXd > __above_lock_fractions
The above lock fractinos for the locked zones and their derivatives.
Definition: BinarySystem.h:136
SEMIMAJOR
SEMIMAJOR
The derivative w.r.t. the semimajor axis in AU.
Definition: DissipationQuantities.h:99
Evolve::BinarySystem::periapsis_evolution_zone_derivs
void periapsis_evolution_zone_derivs(Dissipation::QuantityEntry entry, DissipatingBody &body, unsigned zone_ind, double zone_y_torque_above, double zone_y_torque_below, const std::valarray< Eigen::Vector3d > &zone_torque_deriv, double orbit_y_torque, const std::valarray< Eigen::Vector3d > &orbit_torque_deriv, const std::valarray< Eigen::VectorXd > &above_frac_deriv, double sin_inc, double cos_inc, unsigned locked_zone_ind, double *result) const
Computes the derivatives of the periapsis evolution w.r.t. a zone specific quantity for all zones...
Definition: BinarySystem.cpp:1542
Evolve::DissipatingBody::unlock_zone_spin
void unlock_zone_spin(unsigned zone_index, short direction)
Releases the given zone from a spin-orbit lock.
Definition: DissipatingBody.h:324
Evolve::BinarySystem::inclination_evolution_zone_derivs
void inclination_evolution_zone_derivs(Dissipation::QuantityEntry entry, DissipatingBody &body, unsigned zone_ind, double zone_x_torque_above, double zone_x_torque_below, const std::valarray< Eigen::Vector3d > &zone_torque_deriv, const Eigen::Vector3d &orbit_torque, const std::valarray< Eigen::Vector3d > &orbit_torque_deriv, const std::valarray< Eigen::VectorXd > &above_frac_deriv, double sin_inc, double cos_inc, unsigned locked_zone_ind, double *result) const
Calculates the derivatives of an inclination evolution equation w.r.t. a zone specific quentity...
Definition: BinarySystem.cpp:1462
MOMENT_OF_INERTIA
MOMENT_OF_INERTIA
Definition: DissipationQuantities.h:99
Evolve::DissipatingBody::reached_critical_age
virtual void reached_critical_age(double)
Change the body as necessary at the given age.
Definition: DissipatingBody.h:582
Evolve::BinarySystem::next_stop_age
virtual double next_stop_age() const
The next age when the evolution needs to be stopped for a system change.
Definition: BinarySystem.cpp:2402
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::BinarySystem::above_lock_fraction
double above_lock_fraction(unsigned locked_zone_index, Dissipation::QuantityEntry entry=Dissipation::NO_DERIV, unsigned deriv_zone_index=0, bool secondary_radius=false)
The fraction of an infinitesimal timestep that a zone spends spinning faster than the lock it is in...
Definition: BinarySystem.cpp:2073
Core::Error::Runtime
Any runtime error.
Definition: Error.h:61
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
NUM_DERIVATIVES
NUM_DERIVATIVES
The total number of derivatives supported.
Definition: DissipationQuantities.h:99
Evolve::BinarySystem::__orbit_power
double __orbit_power
The rate at which the orbit gains energy due to tides.
Definition: BinarySystem.h:80
Evolve::BinarySystem::fill_orbit_angmom_gain_deriv
void fill_orbit_angmom_gain_deriv(std::valarray< double > &orbit_angmom_gain_deriv) const
Computes the derivatives w.r.t. the evolution quantities of the orbit angular momentum gain...
Definition: BinarySystem.cpp:1122
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::BinarySystem::reset_evolution
virtual void reset_evolution()
Resets the evolution of the system.
Definition: BinarySystem.cpp:2363
Evolve::DissipatingBody::mass
double mass() const
The mass of the body (constant with age).
Definition: DissipatingBody.h:534
Evolve::BinarySystem::add_body_rate_deriv
void add_body_rate_deriv(const DissipatingBody &body, VALUE_TYPE(DissipatingBody::*func)(Dissipation::QuantityEntry, unsigned, const Eigen::VectorXd &) const, std::valarray< VALUE_TYPE > &orbit_rate_deriv, unsigned offset) const
Adds the derivatives of a rate by which the orbit is changing due to tidal interactions with a single...
Definition: BinarySystem.cpp:1047
Evolve::BinarySystem::age
double age() const
Returns the present age of the system in Gyr.
Definition: BinarySystem.h:833
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::BinarySystem::fill_locked_surface_orbit
void fill_locked_surface_orbit(std::valarray< double > &orbit) const
Implements fill_orbit() for LOCKED_SURFACE_SPIN evolution mode.
Definition: BinarySystem.cpp:1839
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
Evolve::DissipatingBody::spin_jumped
virtual void spin_jumped()
Notifies the body that its spin just discontinously jumped.
Definition: DissipatingBody.h:577
Evolve::BinarySystem::fill_orbit_torque_deriv
void fill_orbit_torque_deriv(Dissipation::QuantityEntry entry, DissipatingBody &body, unsigned zone_ind, std::valarray< Eigen::Vector3d > &orbit_torque_deriv) const
Definition: BinarySystem.cpp:1357
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::BinarySystem::__above_lock_fractions_angmom_deriv
std::valarray< Eigen::VectorXd > __above_lock_fractions_angmom_deriv
The derivatives of the above lock fractions w.r.t. the angular momenta of the zones.
Definition: BinarySystem.h:136
Evolve::BinarySystem::number_locked_zones
unsigned number_locked_zones() const
How many zones on either body are currently locked.
Definition: BinarySystem.h:848
Evolve::BinarySystem::__orbit_power_expansion_error
double __orbit_power_expansion_error
Estimate of the error in __orbit_power due to truncating the tidal potential eccentricity expansion...
Definition: BinarySystem.h:80
Evolve::BinarySystem::__semimajor_rate_evolution
std::list< double > __semimajor_rate_evolution
The rate at which the semimajor axis evolved at each.
Definition: BinarySystem.h:64
Evolve::BinarySystem::__body2
DissipatingBody & __body2
The second body in the system.
Definition: BinarySystem.h:174
Evolve::BinarySystem::reached_critical_age
virtual void reached_critical_age(double age)
Change the system as necessary at the given age.
Definition: BinarySystem.cpp:2396
Evolve::BinarySystem::calculate_above_lock_fractions
void calculate_above_lock_fractions(Eigen::VectorXd &fractions, Dissipation::QuantityEntry entry=Dissipation::NO_DERIV, bool body1_deriv=true)
Calculates the fraction of a timestep above spin-orbit lock for all locked zones. ...
Definition: BinarySystem.cpp:597
Evolve::BinarySystem::__eccentricity
double __eccentricity
The current eccentricity.
Definition: BinarySystem.h:80
Evolve::DissipatingBody::next_stop_age
virtual double next_stop_age() const
The next age when the evolution needs to be stopped for a change in one of the bodies.
Definition: DissipatingBody.h:586
Evolve::BinarySystem::fill_binary_orbit
void fill_binary_orbit(std::valarray< double > &orbit) const
Implements fill_orbit() for BINARY evolution mode.
Definition: BinarySystem.cpp:1855
Evolve::SecondaryDeathCondition
Satisfied when the planet enters below either the roche sphere or the stellar photosphere.
Definition: SecondaryDeathCondition.h:23
Evolve::BinarySystem::__locked_zones
std::list< unsigned > __locked_zones
A list of indices of locked zones.
Definition: BinarySystem.h:129
Evolve::BinarySystem::__orbit_angmom_gain
double __orbit_angmom_gain
The rate at which the orbit gains angular momentum due to tides.
Definition: BinarySystem.h:80
Evolve::BinarySystem::fill_single_orbit
void fill_single_orbit(std::valarray< double > &orbit) const
Implements fill_orbit() for SINGLE evolution mode.
Definition: BinarySystem.cpp:1886
Evolve::BinarySystem::__eccentricity_rate_evolution
std::list< double > __eccentricity_rate_evolution
The rate at which the eccentricity evolved at each.
Definition: BinarySystem.h:64
Evolve::BinarySystem::__evolution_mode
Core::EvolModeType __evolution_mode
The evolution mode from the last call to configure();.
Definition: BinarySystem.h:126
Evolve::BinarySystem::secondary_died
virtual void secondary_died()
Update the system to account for the death of the secondary.
Definition: BinarySystem.cpp:2271
Evolve::BinarySystem::__above_lock_fractions_decomp
Eigen::ColPivHouseholderQR< Eigen::MatrixXd > __above_lock_fractions_decomp
The matrix that defines the problem to solve for __above_lock_fractions.
Definition: BinarySystem.h:166
Evolve::BinarySystem::fill_zone_torque_deriv
void fill_zone_torque_deriv(Dissipation::QuantityEntry entry, DissipatingBody &body, unsigned zone_ind, std::valarray< Eigen::Vector3d > &zone_torque_deriv) const
Calculates the derivatives of the torque on a zone w.r.t. a zone-specific quantity.
Definition: BinarySystem.cpp:1430
Evolve::DissipatingZone::eccentricity_order
virtual unsigned eccentricity_order() const
Definition: DissipatingZone.h:805
Evolve::DissipatingZone::periapsis_evolution
double periapsis_evolution(const Eigen::Vector3d &orbit_torque, const Eigen::Vector3d &zone_torque, Dissipation::QuantityEntry entry=Dissipation::NO_DERIV, const Eigen::Vector3d &orbit_torque_deriv=Eigen::Vector3d(), const Eigen::Vector3d &zone_torque_deriv=Eigen::Vector3d())
The rate at which the periapsis of the orbit/reference zone in this zone&#39;s equatorial plane is changi...
Definition: DissipatingZone.cpp:658
ECCENTRICITY
ECCENTRICITY
The derivative w.r.t. the eccentricity.
Definition: DissipationQuantities.h:88
Evolve::BinarySystem::spin_angmom_evolution_zone_derivs
void spin_angmom_evolution_zone_derivs(Dissipation::QuantityEntry entry, DissipatingBody &body, unsigned zone_ind, double zone_z_torque_above, double zone_z_torque_below, const std::valarray< Eigen::Vector3d > &zone_torque_deriv, const std::valarray< Eigen::VectorXd > &above_frac_deriv, unsigned locked_zone_ind, double *result) const
Definition: BinarySystem.cpp:1635
Evolve::BinarySystem::jacobian
int jacobian(double age, const double *parameters, Core::EvolModeType evolution_mode, double *param_derivs, double *age_derivs)
Definition: BinarySystem.cpp:2146
Evolve::DissipatingZone::set_evolution_rates
void set_evolution_rates(double angular_momentum, double inclination, double periapsis)
Set evolution rates for this zone&#39;s properties for storing in the eveloution.
Definition: DissipatingZone.h:360
Evolve::DissipatingZone::set_reference_zone_angmom
void set_reference_zone_angmom(double reference_angmom)
Defines the angular momentum of the reference zone for single body evolution.
Definition: DissipatingZone.h:395
Evolve::BinarySystem::single_body_jacobian
void single_body_jacobian(double *param_derivs, double *age_derivs) const
Jacobian for the evolution of the rotation of the zones of body 0 if no other body is present...
Definition: BinarySystem.cpp:394
Evolve::BinarySystem::__above_lock_fractions_inclination_deriv
std::valarray< Eigen::VectorXd > __above_lock_fractions_inclination_deriv
The derivatives of the above lock fractions w.r.t. the inclinations of the zones. ...
Definition: BinarySystem.h:136
Evolve::BinarySystem::binary_jacobian
void binary_jacobian(double *param_derivs, double *age_derivs) const
Calculates the jacobian for the evolution of a system of two bodies.
Definition: BinarySystem.cpp:1679
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::BinarySystem::fill_orbit_power_deriv
void fill_orbit_power_deriv(std::valarray< double > &orbit_power_deriv) const
Computes the derivatives w.r.t. the evolution quantities of the orbit energy gain.
Definition: BinarySystem.cpp:1107
BinarySystem.h
Defines the BinarySystem class.
Evolve::BinarySystem::rewind_evolution
virtual void rewind_evolution(unsigned nsteps)
Discards the last steps from the evolution.
Definition: BinarySystem.cpp:2373
Evolve::BinarySystem::evolution_mode
Core::EvolModeType evolution_mode()
The evolution mode of last call to configure().
Definition: BinarySystem.h:996
Evolve::BinarySystem::update_above_lock_fractions
void update_above_lock_fractions()
Solves for and sets the above lock fractions and their derivatives.
Definition: BinarySystem.cpp:821
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
Evolve::ZoneOrientation::inclination
double inclination(bool evolution_rate=false) const
The angle between the angular momenta of the zone and the orbit.
Definition: ZoneOrientation.h:54
NO_DERIV
NO_DERIV
The quantity itself, undifferentiated.
Definition: DissipationQuantities.h:63
Evolve::BinarySystem::__orbit_torque_expansion_error
Eigen::Vector3d __orbit_torque_expansion_error
An estiamte of the error in ::__orbit_torque due to truncating the eccentricity series of the tidal p...
Definition: BinarySystem.h:119
Evolve::BinarySystem::eccentricity_jacobian
void eccentricity_jacobian(const std::valarray< double > &orbit_power_deriv, const std::valarray< double > &orbit_angmom_gain_deriv, bool a6p5, double *param_derivs, double &age_deriv) const
Computes the row of the jacobian corresponding to the eccentricity.
Definition: BinarySystem.cpp:1182
Core::EvolModeType
EvolModeType
The various evolution modes.
Definition: Common.h:42
ORBITAL_FREQUENCY
ORBITAL_FREQUENCY
The derivative w.r.t. the orbital frequency.
Definition: DissipationQuantities.h:63
Evolve::BinarySystem::__above_lock_fractions_inertia_deriv
std::valarray< Eigen::VectorXd > __above_lock_fractions_inertia_deriv
The derivatives of the above lock fractions w.r.t. the moments of inertia of the zones.
Definition: BinarySystem.h:136
Evolve::DissipatingBody::lock_zone_spin
void lock_zone_spin(unsigned zone_index, int orbital_frequency_multiplier, int spin_frequency_multiplier)
Definition: DissipatingBody.h:314
Evolve::BinarySystem::semimajor
double semimajor() const
The current semimajor axis of the system.
Definition: BinarySystem.h:852
SPIN_ANGMOM
SPIN_ANGMOM
The derivative w.r.t. the spin angular momentum in .
Definition: DissipationQuantities.h:99
Evolve::DissipatingZone::inclination_evolution
double inclination_evolution(const Eigen::Vector3d &orbit_torque, const Eigen::Vector3d &zone_torque, Dissipation::QuantityEntry entry=Dissipation::NO_DERIV, const Eigen::Vector3d &orbit_torque_deriv=Eigen::Vector3d(), const Eigen::Vector3d &zone_torque_deriv=Eigen::Vector3d())
The rate at which the inclination between this zone and the orbit is changing.
Definition: DissipatingZone.cpp:762
Evolve::BinarySystem::eccentricity_evolution_expansion_error
double eccentricity_evolution_expansion_error() const
Estimate of the error in the value returned by eccentricity_evolution() due to truncating the tidal p...
Definition: BinarySystem.cpp:485
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::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::BinarySystem::locked_surface_jacobian
void locked_surface_jacobian(double *param_derivs, double *age_derivs) const
Jacobian for the evolution of the rotation of the zones of body 0 with the topmost zone rotating with...
Definition: BinarySystem.cpp:68
Evolve::BinarySystem::__body1
DissipatingBody & __body1
The first body in the system.
Definition: BinarySystem.h:174
Evolve::BinarySystem::above_lock_fractions_deriv
Eigen::VectorXd above_lock_fractions_deriv(Dissipation::QuantityEntry entry, DissipatingBody &body, unsigned zone_index)
Calculates derivatives of the above lock fractions w.r.t. quantities of non-surface zones...
Definition: BinarySystem.cpp:700
Evolve::BinarySystem::__above_lock_fractions_periapsis_deriv
std::valarray< Eigen::VectorXd > __above_lock_fractions_periapsis_deriv
The derivatives of the above lock fractions w.r.t. the periapses of the zones.
Definition: BinarySystem.h:136
Evolve::CombinedStoppingCondition
A class combining the the outputs of multiple stopping conditions.
Definition: CombinedStoppingCondition.h:21
Evolve::BinarySystem::__eccentricity_rate
double __eccentricity_rate
The current rate at which the eccentricity is changing.
Definition: BinarySystem.h:111
Evolve::DissipatingZone::dissipative
virtual bool dissipative() const =0
Should return true iff the zone has some non-zero dissipation.
Evolve::BinarySystem::number_zones
unsigned number_zones() const
The total number of zones in both system bodies.
Definition: BinarySystem.h:842
Evolve::BinarySystem::above_lock_problem_deriv_correction
void above_lock_problem_deriv_correction(Dissipation::QuantityEntry entry, bool body1_deriv, Eigen::MatrixXd &matrix, Eigen::VectorXd &rhs) const
Makes corrections to the matrix and RHS to accomodate the given derivative for the linear problem tha...
Definition: BinarySystem.cpp:501
Evolve::BinarySystem::semimajor_evolution_expansion_error
double semimajor_evolution_expansion_error() const
Estimate of the error in the value returned by semimajor_evolution() due to truncating the tidal pote...
Definition: BinarySystem.h:294
Evolve::BinarySystem::angle_evolution_orbit_deriv
void angle_evolution_orbit_deriv(Dissipation::QuantityEntry entry, double angmom_deriv, DissipatingBody &body, unsigned zone_ind, double sin_inc, double cos_inc, unsigned locked_zone_ind, double &inclination, double &periapsis) const
Computes the partial derivatives w.r.t. semimamjor axis or eccentricity of the inclination and periap...
Definition: BinarySystem.cpp:1287
Evolve::BinarySystem::semimajor_jacobian
void semimajor_jacobian(const std::valarray< double > &orbit_power_deriv, bool a6p5, double *param_derivs, double &age_deriv) const
Computes the row of the jacobian corresponding to the semimajor axis.
Definition: BinarySystem.cpp:1165
Evolve::BinarySystem::release_lock
virtual void release_lock(unsigned locked_zone_index, short direction)
Releases the lock to one of the locked zones.
Definition: BinarySystem.cpp:2321
Evolve::BinarySystem::__orbital_angmom
double __orbital_angmom
The current orbital angular momentum.
Definition: BinarySystem.h:80
Evolve::BinarySystem::eccentricity
double eccentricity() const
The current eccentricity of the system.
Definition: BinarySystem.h:855
Evolve::BinarySystem::stopping_conditions
virtual CombinedStoppingCondition * stopping_conditions()
Conditions detecting the next possible doscontinuity in the evolution.
Definition: BinarySystem.cpp:2385
Evolve::BinarySystem::configure
virtual int configure(bool initialize, double age, double semimajor, double eccentricity, const double *spin_angmom, const double *inclination, const double *periapsis, Core::EvolModeType evolution_mode)
Sets the current state of the system.
Definition: BinarySystem.cpp:1912
Evolve::BinarySystem::__orbital_energy
double __orbital_energy
The current orbital energy.
Definition: BinarySystem.h:80
Evolve::BinarySystem::fill_orbit_torque_and_power
void fill_orbit_torque_and_power()
Update the values of ::__orbit_power, ::__orbit_torque, ::__orbit_angmom_gain and their associated ex...
Definition: BinarySystem.cpp:849
Evolve::BinarySystem::eccentricity_order
virtual unsigned eccentricity_order() const
Definition: BinarySystem.cpp:2408
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::BinarySystem::locked_surface_differential_equations
int locked_surface_differential_equations(double *evolution_rates, bool expansion_error) const
Differential equations for the rotation of the zones of body 0 with the topmost zone rotating with a ...
Definition: BinarySystem.cpp:31
Evolve::BinarySystem::fill_orbit
Core::EvolModeType fill_orbit(std::valarray< double > &orbit) const
Fills an array with the parameters expected by differential_equations() and jacobian(), returning the evolution mode.
Definition: BinarySystem.cpp:2062
Evolve::DissipatingZone::spin_frequency
double spin_frequency() const
The spin frequency of the given zone.
Definition: DissipatingZone.h:584
Evolve::DissipatingBody::number_zones
virtual unsigned number_zones() const =0
The number of zones the body consists of.
Evolve::BinarySystem::angle_evolution_age_deriv
void angle_evolution_age_deriv(DissipatingBody &body, unsigned zone_ind, double sin_inc, double cos_inc, unsigned locked_zone_ind, double &inclination, double &periapsis) const
Computes the partial age derivatives of the inclination or periapsis evolutiono equations.
Definition: BinarySystem.cpp:1209