Planetary Orbital Evolution due to Tides
Orbital evolution of two objects experiencing tides
Common.cpp
Go to the documentation of this file.
1 
8 #include "Common.h"
9 
10 namespace std {
11  std::ostream &operator<<(std::ostream &os,
12  const std::valarray<double> &arr)
13  {
14  for(size_t i=0; i<arr.size(); i++) {
15  if(i) os << ", ";
16  os << arr[i];
17  }
18  return os;
19  }
20 
21  std::ostream &operator<<(std::ostream &os, const std::list<double> &l)
22  {
23  for(
24  std::list<double>::const_iterator i=l.begin();
25  i!=l.end();
26  ++i
27  )
28  os << *i << ", ";
29  return os;
30  }
31 }
32 
33 namespace Core {
34 
35  std::ostream &operator<<(std::ostream &os,
36  const Core::EvolModeType &evol_mode)
37  {
38  switch(evol_mode) {
39  case BINARY: os << "BINARY"; break;
40  case SINGLE: os << "SINGLE"; break;
41  case LOCKED_SURFACE_SPIN: os << "LOCKED_SURFACE_SPIN"; break;
42  case TABULATION: os << "TABULATION";
43  }
44  return os;
45  }
46 
47  const double NUMERIC_SAFETY=100.0*std::numeric_limits<double>::epsilon();
48 
49  std::valarray<double> solve_cubic(double c0,
50  double c1,
51  double c2,
52  double c3) {
53  double x[3];
54  int n = 0;
55  if (c3 == 0)
56  n = gsl_poly_solve_quadratic(c2, c1, c0, &x[0], &x[1]);
57  else
58  n = gsl_poly_solve_cubic(c2/c3,
59  c1/c3,
60  c0/c3,
61  &x[0],
62  &x[1],
63  &x[2]);
64  std::valarray<double> result(x, n);
65  return result;
66  }
67 
68  double estimate_zerocrossing(double x0, double y0,
69  double x1, double y1,
70  double dy0, double dy1)
71  {
72  double dx=x1-x0, dy=y1-y0, slope_inv=dx/dy;
73  double linear_solution=x0 - y0*slope_inv;
74  if(linear_solution<x0 || linear_solution>x1)
75  linear_solution=x1 - y1*slope_inv;
76  if(std::isnan(dy0) || std::isnan(dy1)) {
77 #ifndef NDEBUG
78  /* std::cerr << "Linear zerocrossing between ("
79  << x0 << ", " << y0
80  << ") and (" << x1 << ", " << y1 << ")="
81  << linear_solution
82  << std::endl;*/
83 #endif
84  return linear_solution;
85  } else {
86  double x1_2=std::pow(x1, 2), x0_2=std::pow(x0, 2),
87  a=(dx*(dy1+dy0) - 2.0*dy)/std::pow(dx, 3),
88  b=(dy1-dy0 - 3.0*a*(x1_2-x0_2))/(2.0*dx),
89  c=dy0 - 3.0*a*x0_2 - 2.0*b*x0,
90  d=y0 - a*x0*x0_2 - b*x0_2 - c*x0;
91  std::valarray<double> solutions=solve_cubic(d, c, b, a);
92 #ifndef NDEBUG
93  /* std::cerr << "Cubic zerocrossing between ("
94  << x0 << ", " << y0
95  << ", " << dy0 << ") and ("
96  << x1 << ", " << y1 << ", " << dy1
97  << "), coef=(" << a << ", " << b << ", " << c << ", "
98  << d << "), solutions=(" << solutions << ")";*/
99 #endif
100  for(size_t i=0; i<solutions.size(); i++)
101  if(solutions[i]>=x0 && solutions[i]<=x1) {
102 #ifndef NDEBUG
103  /* std::cerr << ", selected: "
104  << solutions[i] << std::endl;*/
105 #endif
106  return solutions[i];
107  }
108  if(y0*y1<=0) {
109 #ifndef NDEBUG
110  /* std::cerr << ", fallback to linear: "
111  << linear_solution
112  << std::endl;*/
113 #endif
114  return linear_solution;
115  }
116  std::ostringstream msg;
117  msg.precision(16);
118  msg << "No solution to "
119  << a << "*x**3 + "
120  << b << "*x**2 + "
121  << c << "*x + "
122  << d
123  << " (x0=" << x0
124  << ", y0=" << y0
125  << ", x1=" << x1
126  << ", y1=" << y1
127  << ", dy0=" << dy0
128  << ", dy1=" << dy1
129  << " found in the given range in estimate zerocrossing. "
130  "Solutions(" << solutions.size() << "): ";
131  if(solutions.size()>0) msg << solutions[0];
132  for(size_t i=1; i<solutions.size(); i++)
133  msg << ", " << solutions[i];
134  throw Error::BadFunctionArguments(msg.str());
135  }
136  }
137 
138  double quadratic_zerocrossing(double x0, double y0,
139  double x1, double y1,
140  double x2, double y2,
141  double require_range_low,
142  double require_range_high)
143  {
144  double d12=(y1-y2)/(x1-x2), d02=(y0-y2)/(x0-x2), dx10=(x1-x0),
145  a=d12/(x1-x0) - d02/dx10,
146  b=d12 + d02*(x1+x2)/dx10 - d12*(x1+x2)/dx10,
147  c=y2 + d12*x1*x2/dx10 - d12*x2 - d02*x1*x2/dx10,
148  xpre, ypre, xpost, ypost;
149  if(y0*y1<0) {xpre=x0; xpost=x1; ypre=y0; ypost=y1;}
150  else if(y1*y2<0) {xpre=x1; xpost=x2; ypre=y1; ypost=y2;}
151  else {
152  std::ostringstream msg;
153  msg.precision(16);
154  msg << "No sign change found in quadratic_zerocrossing(x0="
155  << x0 << ", y0=" << y0 << ", x1=" << x1 << ", y1=" << y1
156  << ", x2=" << x2 << ", y2=" << y2 << ")";
157  throw Error::BadFunctionArguments(msg.str());
158  }
159  std::valarray<double> solutions=solve_cubic(c, b, a, 0);
160  if(std::isnan(require_range_low)) {
161  assert(std::isnan(require_range_high));
162  require_range_low=x0;
163  require_range_high=x2;
164  }
165 #ifndef NDEBUG
166  /* std::cerr << "Quadratic zerocrossing between (" << x0 << ", "
167  << y0
168  << "), (" << x1 << ", " << y1 << "), (" << x2 << ", " << y2
169  << ") in range (" << require_range_low << ", "
170  << require_range_high
171  << "), coef=(" << a << ", " << b << ", " << c << "), solutions=("
172  << solutions << ")";*/
173 #endif
174  for(size_t i=0; i<solutions.size(); i++)
175  if(solutions[i]>=require_range_low &&
176  solutions[i]<=require_range_high) {
177 #ifndef NDEBUG
178  /* std::cerr << ", selected: " << solutions[i]
179  << std::endl;*/
180 #endif
181  return solutions[i];
182  }
183 #ifndef NDEBUG
184  // std::cerr << ", fallback to linear: ";
185 #endif
186  return estimate_zerocrossing(xpre, ypre, xpost, ypost);
187  }
188 
189  double cubic_zerocrossing(double x0, double y0,
190  double x1, double y1,
191  double x2, double y2,
192  double x3, double y3,
193  double require_range_low,
194  double require_range_high)
195  {
196  double xpre, ypre, xpost, ypost;
197  if(y0*y1<=0) {xpre=x0; xpost=x1; ypre=y0; ypost=y1;}
198  else if(y1*y2<=0) {xpre=x1; xpost=x2; ypre=y1; ypost=y2;}
199  else if(y2*y3<=0) {xpre=x2; xpost=x3; ypre=y2; ypost=y3;}
200  else {
201  std::ostringstream msg;
202  msg.precision(16);
203  msg << "No sign change found in cubic_zerocrossing(x0="
204  << x0 << ", y0=" << y0 << ", x1=" << x1 << ", y1=" << y1
205  << ", x2=" << x2 << ", y2=" << y2
206  << ", x3=" << x3 << ", y3=" << y3 << ")";
207  throw Error::BadFunctionArguments(msg.str());
208  }
209  gsl_vector *cubic_coef=cubic_coefficients(x0, y0,
210  x1, y1,
211  x2, y2,
212  x3, y3);
213  std::valarray<double> solutions=solve_cubic(
214  gsl_vector_get(cubic_coef, 0),
215  gsl_vector_get(cubic_coef, 1),
216  gsl_vector_get(cubic_coef, 2),
217  gsl_vector_get(cubic_coef, 3));
218  if(std::isnan(require_range_low)) {
219  assert(std::isnan(require_range_high));
220  require_range_low=x0;
221  require_range_high=x3;
222  }
223 #ifndef NDEBUG
224  /* std::cerr << "Cubic zerocrossing between (" << x0 << ", " << y0
225  << "), (" << x1 << ", " << y1 << "), (" << x2 << ", " << y2
226  << "), (" << x3 << ", " << y3 << ") in range ("
227  << require_range_low << ", " << require_range_high
228  << "), coef=(";
229  for(int i=3; i>=0; i--) {
230  std::cerr << gsl_vector_get(cubic_coef, i);
231  if(i>0) std::cerr << ", ";
232  }
233  std::cerr << "), solutions=(" << solutions << ")";*/
234 #endif
235  gsl_vector_free(cubic_coef);
236  for(size_t i=0; i<solutions.size(); i++)
237  if(solutions[i]>=require_range_low &&
238  solutions[i]<=require_range_high) {
239 #ifndef NDEBUG
240  /* std::cerr << ", selected: " << solutions[i]
241  << std::endl;*/
242 #endif
243  return solutions[i];
244  }
245 #ifndef NDEBUG
246  // std::cerr << ", fallback to linear: ";
247 #endif
248  return estimate_zerocrossing(xpre, ypre, xpost, ypost);
249  }
250 
251  double estimate_extremum(double x0, double y0,
252  double x1, double y1,
253  double dy0, double dy1,
254  double *extremum_y)
255  {
256  assert(dy0*dy1<0);
257  double dx=x1-x0;
258  double x1_2=std::pow(x1, 2), x0_2=std::pow(x0, 2),
259  a=(dx*(dy1+dy0) - 2.0*(y1-y0))/std::pow(dx, 3),
260  b=(dy1-dy0 - 3.0*a*(x1_2-x0_2))/(2.0*dx),
261  c=dy0 - 3.0*a*x0_2 - 2.0*b*x0,
262  d=y0 - a*x0*x0_2 - b*x0_2 - c*x0,
263  linear_extremum_x=x0 - dy0*dx/(dy1-dy0),
264  linear_extremum_y=(linear_extremum_x-x0<x1-linear_extremum_x ?
265  y0+dy0*(linear_extremum_x-x0) :
266  y1-dy1*(x1-linear_extremum_x));
267  std::valarray<double> solutions=solve_cubic(c, 2.0*b, 3.0*a, 0);
268 #ifndef NDEBUG
269  /* std::cerr << "Cubic extrema between (" << x0 << ", "
270  << y0 << ", " << dy0
271  << ") and (" << x1 << ", " << y1 << ", " << dy1 << "), coef=("
272  << a << ", " << b << ", " << c << ", " << d
273  << "), solutions(" << solutions.size() << "): (";
274  if(solutions.size()==0) std::cerr << "no solutions)";
275  else {
276  std::cerr << solutions[0] << ", "
277  << a*std::pow(solutions[0], 3) + b*std::pow(solutions[0], 2) +
278  c*solutions[0] + d << "), (" << solutions[1] << ", "
279  << a*std::pow(solutions[1], 3) + b*std::pow(solutions[1], 2) +
280  c*solutions[1] + d << ")";
281  }*/
282 #endif
283  for(size_t i=0; i<solutions.size(); i++) {
284  if(solutions[i]>x0 && solutions[i]<x1) {
285  double x=solutions[i];
286  if(extremum_y!=NULL) {
287  double ax3=a*x*x*x, bx2=b*x*x, cx=c*x;
288  *extremum_y=ax3 + bx2 + cx + d;
289  if(std::abs(*extremum_y)/std::max(
290  std::max(std::abs(ax3), std::abs(bx2)),
291  std::max(std::abs(cx), std::abs(d)))<1e-10)
292  *extremum_y=linear_extremum_y;
293  }
294 #ifndef NDEBUG
295  /* std::cerr << ", selected: (" << x;
296  if(extremum_y) std::cerr << ", " << *extremum_y;
297  std::cerr << ")" << std::endl;*/
298 #endif
299  return x;
300  }
301  }
302  if(extremum_y!=NULL) *extremum_y=linear_extremum_y;
303 #ifndef NDEBUG
304  /* std::cerr << ", falling back to linear: ("
305  << linear_extremum_x << ", "
306  << linear_extremum_y << ")" << std::endl;*/
307 #endif
308  return linear_extremum_x;
309  }
310 
312  bool indistinguishable(double x0, double y0, double x1, double y1)
313  {
314  return (std::abs(x1-x0)/std::max(std::abs(x1), std::abs(x0))
315  <NUMERIC_SAFETY
316  &&
317  std::abs(y1-y0)/std::max(std::abs(y1), std::abs(y0))
318  <NUMERIC_SAFETY);
319  }
320 
321  double quadratic_extremum(double x0, double y0, double x1,
322  double y1, double x2, double y2, double *extremum_y)
323  {
324 #ifndef NDEBUG
325  assert((y1-y0)*(y2-y1)<=0);
326 #endif
327  if(indistinguishable(x0, y0, x1, y1) ||
328  indistinguishable(x1, y1, x2, y2)) {
329 #ifndef NDEBUG
330  /* std::cerr << "Quandratic extremum with indistinguisable two points. "
331  "Returing middle point as answer." << std::endl;*/
332 #endif
333  if(extremum_y) *extremum_y=y1;
334  return x1;
335  }
336  double s02=(y0-y2)/(x0-x2), s12=(y1-y2)/(x1-x2), a=(s02 - s12)/(x0-x1),
337  extremum_x=0.5*(x0 + x2 - s02/a);
338  if(extremum_y)
339  *extremum_y=y0 - s02*s02/(2.0*a) + s02*x2
340  - a*(x0*x0 + x2*x2)/2.0;
341 #ifndef NDEBUG
342  /* double b=-2.0*a*extremum_x;
343  std::cerr << "Quadratic extremum between (" << x0 << ", "
344  << y0 << "), ("
345  << x1 << ", " << y1 << "), (" << x2 << ", " << y2 << "), coef=("
346  << a << ", " << b << ", "
347  << *extremum_y - a*extremum_x*extremum_x - b*extremum_x
348  << "): (" << extremum_x << ", " << *extremum_y << ")"
349  << std::endl;*/
350 #endif
351  return extremum_x;
352  }
353 
354  double cubic_extremum(double x0, double y0, double x1,
355  double y1, double x2,double y2, double x3, double y3,
356  double *extremum_y, double require_range_low,
357  double require_range_high)
358  {
359 #ifndef NDEBUG
360  assert((y1-y0)*(y2-y1)<=0 || (y2-y1)*(y3-y2)<=0);
361 #endif
362  double extremum_x;
363  if(indistinguishable(x0, y0, x1, y1)) {
364 #ifndef NDEBUG
365  /* std::cerr << "Cubic extremum with indistinguishable first two "
366  "points." << std::endl;*/
367 #endif
368  if((y2-y1)*(y3-y2)<=0)
369  extremum_x=quadratic_extremum(x1,y1, x2,y2, x3,y3, extremum_y);
370  else {
371  extremum_x=x1;
372  if(extremum_y) *extremum_y=y1;
373  }
374  } else if(indistinguishable(x1, y1, x2, y2)) {
375 #ifndef NDEBUG
376  /* std::cerr << "Cubic extremum with indistinguishable middle two "
377  "points." << std::endl;*/
378 #endif
379  double xmid=0.5*(x1+x2), ymid=0.5*(y1+y2);
380  if((ymid-y0)*(y3-ymid)<=0)
381  extremum_x=quadratic_extremum(x0, y0, xmid, ymid, x3, y3,
382  extremum_y);
383  else if((y1-y0)*(y2-y1)<=0) {
384  extremum_x=x1;
385  if(extremum_y) *extremum_y=y1;
386  } else {
387  extremum_x=x2;
388  if(extremum_y) *extremum_y=y2;
389  }
390  } else if(indistinguishable(x2, y2, x3, y3)) {
391 #ifndef NDEBUG
392  /* std::cerr << "Cubic extremum with indistinguishable last two "
393  "points." << std::endl;*/
394 #endif
395  if((y1-y0)*(y2-y1)<=0)
396  extremum_x=quadratic_extremum(x0,y0, x1,y1, x2,y2, extremum_y);
397  else {
398  extremum_x=x2;
399  if(extremum_y) *extremum_y=y2;
400  }
401  } else {
402  gsl_vector *cubic_coef=
403  cubic_coefficients(x0,y0, x1,y1, x2,y2, x3,y3);
404  double a=3.0*gsl_vector_get(cubic_coef, 3),
405  b=2.0*gsl_vector_get(cubic_coef, 2),
406  c=gsl_vector_get(cubic_coef, 1), sqrtD=std::sqrt(b*b-4.0*a*c);
407 #ifndef NDEBUG
408  /* double extremum_x1=(-b-sqrtD)/(2.0*a),
409  extremum_x2=(-b+sqrtD)/(2.0*a),
410  d=y0 - a*x0*x0*x0/3.0 - b*x0*x0/2.0 - c*x0;*/
411  if(std::isnan(require_range_low)) {
412  assert(std::isnan(require_range_high));
413  require_range_low=x0;
414  require_range_high=x3;
415  }
416  /* std::cerr << "Cubic extrema between (" << x0 << ", "
417  << y0 << "), ("
418  << x1 << ", " << y1 << "), (" << x2 << ", " << y2
419  << "), ("
420  << x3 << ", " << y3 << ") in range ("
421  << require_range_low
422  << ", " << require_range_high << "), coef=("
423  << a/3.0 << ", " << b/2.0 << ", " << c << ", " << d
424  << "), solutions: (" << extremum_x1 << ", "
425  << (a*std::pow(extremum_x1,3)/3.0
426  +
427  b*std::pow(extremum_x1,2)/2.0
428  +
429  c*extremum_x1 + d)
430  << ") and (" << extremum_x2 << ", "
431  << a*std::pow(extremum_x2,3)/3.0 +
432  b*std::pow(extremum_x2,2)/2.0 + c*extremum_x2 + d
433  << ")";*/
434 #endif
435  if(a<0) {a*=-1; b*=-1; c*=-1;}
436  extremum_x=(-b-sqrtD)/(2.0*a);
437  if(extremum_x<require_range_low || extremum_x>require_range_high)
438  extremum_x=(-b+sqrtD)/(2.0*a);
439  if(extremum_y) {
440  *extremum_y=0;
441  double xpow=1.0;
442  for(size_t i=0; i<4; i++) {
443  *extremum_y+=xpow*gsl_vector_get(cubic_coef, i);
444  xpow*=extremum_x;
445  }
446  }
447  gsl_vector_free(cubic_coef);
448  }
449  if(extremum_x<require_range_low || extremum_x>require_range_high) {
450 #ifndef NDEBUG
451  // std::cerr << ", fallback to quadratic: ";
452 #endif
453  if((y1-y0)*(y2-y1)<=0)
454  extremum_x=quadratic_extremum(x0, y0,
455  x1, y1,
456  x2, y2,
457  extremum_y);
458  if(
459  (
460  extremum_x<require_range_low
461  ||
462  extremum_x>require_range_high
463  )
464  &&
465  (y2 - y1) * (y3 - y2) <= 0
466  )
467  extremum_x = quadratic_extremum(x1, y1,
468  x2, y2,
469  x3, y3,
470  extremum_y);
471  if(extremum_x<require_range_low || extremum_x>require_range_high) {
472  std::ostringstream msg;
473  msg.precision(16);
474  msg.setf(std::ios_base::scientific);
475  msg << "No extremum found in the range ("
476  << require_range_low
477  << ", " << require_range_high << ") from points ("
478  << x0 << ", " << y0 << "), ("
479  << x1 << ", " << y1 << "), ("
480  << x2 << ", " << y2 << ") and ("
481  << x3 << ", " << y3 << ")!";
482  throw Error::BadFunctionArguments(msg.str());
483  }
484  }
485 #ifndef NDEBUG
486  /* std::cerr << ", selected: (" << extremum_x;
487  if(extremum_y) std::cerr << ", " << *extremum_y;
488  std::cerr << ")" << std::endl;*/
489 #endif
490  return extremum_x;
491  }
492 
493  gsl_vector *cubic_coefficients(double x0, double y0,
494  double x1, double y1,
495  double x2, double y2,
496  double x3, double y3)
497  {
498  std::list<double> xs, ys;
499  xs.push_back(x0); ys.push_back(y0);
500  xs.push_back(x1); ys.push_back(y1);
501  xs.push_back(x2); ys.push_back(y2);
502  xs.push_back(x3); ys.push_back(y3);
503  return polynomial_coefficients(xs.begin(), ys.begin(),4);
504  }
505 
506 } //End Core namespace.
Core::Error::BadFunctionArguments
Function arguments do not satisfy some requirement.
Definition: Error.h:73
Core::cubic_coefficients
gsl_vector * cubic_coefficients(double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3)
Returns the polynomial coefficients of the cubic going through four points.
Definition: Common.cpp:493
Core::estimate_zerocrossing
double estimate_zerocrossing(double x0, double y0, double x1, double y1, double dy0, double dy1)
Finds a zero of a smooth curve defined by the given points and derivatives.
Definition: Common.cpp:68
std
Definition: Common.cpp:10
Core::solve_cubic
std::valarray< double > solve_cubic(double c0, double c1, double c2, double c3)
Solves the cubic equation .
Definition: Common.cpp:49
Core::cubic_zerocrossing
double cubic_zerocrossing(double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3, double require_range_low, double require_range_high)
Finds the abscissa of a zero of a cubic defined by four points.
Definition: Common.cpp:189
std::operator<<
std::ostream & operator<<(std::ostream &os, const std::valarray< double > &arr)
Outputs a valarray as a sequence of &#39;, &#39; separated values.
Definition: Common.cpp:11
Core::indistinguishable
bool indistinguishable(double x0, double y0, double x1, double y1)
Are the two points indistinguishable (up to NUMERIC_SAFETY).
Definition: Common.cpp:312
Core::estimate_extremum
double estimate_extremum(double x0, double y0, double x1, double y1, double dy0, double dy1, double *extremum_y)
Finds an extremum of a cubic defined by two points and derivatives.
Definition: Common.cpp:251
Core::LOCKED_SURFACE_SPIN
The surface rotation of one of the bodies is locked to a prescribed value.
Definition: Common.h:45
Core::polynomial_coefficients
gsl_vector * polynomial_coefficients(ITERATOR x_i, ITERATOR y_i, size_t num_points)
Finds the polynomial that goes through the points iterated over by x_i and y_i.
Definition: Common.h:233
Core::SINGLE
There is only one body in the system (only its rotation evolves).
Definition: Common.h:51
Core::quadratic_extremum
double quadratic_extremum(double x0, double y0, double x1, double y1, double x2, double y2, double *extremum_y)
Finds the abscissa of the extremum of the quadratic function passing through three points...
Definition: Common.cpp:321
Core::cubic_extremum
double cubic_extremum(double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3, double *extremum_y, double require_range_low, double require_range_high)
Finds the abscissa of an extremum of the cubic function passing through four points.
Definition: Common.cpp:354
Common.h
Declaration of some general purpose utilities.
Core::EvolModeType
EvolModeType
The various evolution modes.
Definition: Common.h:42
Core::BINARY
Two bodies orbiting each other.
Definition: Common.h:48
Core::quadratic_zerocrossing
double quadratic_zerocrossing(double x0, double y0, double x1, double y1, double x2, double y2, double require_range_low, double require_range_high)
Finds the abscissa of a zero of a quadratic defined by three points.
Definition: Common.cpp:138