{"action":"created","comment":{"url":"https://api.github.com/repos/mockingbirdnest/Principia/pulls/comments/30839415","id":30839415,"diff_hunk":"@@ -0,0 +1,209 @@\n+﻿#pragma once\n+\n+#include \"integrators/explicit_embedded_runge_kutta_nyström_integrator.hpp\"\n+\n+#include <algorithm>\n+#include <cmath>\n+#include <ctime>\n+#include <vector>\n+\n+#include \"geometry/sign.hpp\"\n+#include \"glog/logging.h\"\n+#include \"quantities/quantities.hpp\"\n+\n+namespace principia {\n+\n+using geometry::Sign;\n+using quantities::Difference;\n+using quantities::Quotient;\n+\n+namespace integrators {\n+\n+inline ExplicitEmbeddedRungeKuttaNyströmIntegrator const&\n+DormandElMikkawyPrince1986RKN434FM() {\n+  static ExplicitEmbeddedRungeKuttaNyströmIntegrator const integrator(\n+      // c\n+      {0.0, 1.0 / 4.0, 7.0 / 10.0, 1.0},\n+      // a\n+      {{},\n+       {1.0 /   32.0},\n+       {7.0 / 1000.0, 119.0 / 500.0},\n+       {1.0 /   14.0,   8.0 /  27.0, 25.0 / 189.0}},\n+      // b̂\n+      { 1.0 /  14.0,   8.0 /  27.0,  25.0 / 189.0,  0.0},\n+      // b̂′\n+      { 1.0 /  14.0,  32.0 /  81.0, 250.0 / 567.0,  5.0 / 54.0},\n+      // b\n+      {-7.0 / 150.0,  67.0 / 150.0,   3.0 /  20.0, -1.0 / 20.0},\n+      // b′\n+      {13.0 /  21.0, -20.0 /  27.0, 275.0 / 189.0, -1.0 /  3.0},\n+      3);  // lower_order\n+  return integrator;\n+}\n+\n+inline ExplicitEmbeddedRungeKuttaNyströmIntegrator::\n+           ExplicitEmbeddedRungeKuttaNyströmIntegrator(\n+    std::vector<double> const& c,\n+    std::vector<std::vector<double>> const& a,\n+    std::vector<double> const& b_hat,\n+    std::vector<double> const& b_prime_hat,\n+    std::vector<double> const& b,\n+    std::vector<double> const& b_prime,\n+    int const lower_order)\n+    : stages_(c.size()),\n+      lower_order_(lower_order),\n+      c_(c),\n+      a_(a),\n+      b_hat_(b_hat),\n+      b_prime_hat_(b_prime_hat),\n+      b_(b),\n+      b_prime_(b_prime) {\n+  CHECK_EQ(a_.size(), stages_);\n+  for(int i = 0; i < stages_; ++i) {\n+    CHECK_EQ(a_[i].size(), i);\n+  }\n+  CHECK_EQ(b_hat_.size(), stages_);\n+  CHECK_EQ(b_prime_hat_.size(), stages_);\n+  CHECK_EQ(b_.size(), stages_);\n+  CHECK_EQ(b_prime_.size(), stages_);\n+}\n+\n+template<typename Position>\n+void ExplicitEmbeddedRungeKuttaNyströmIntegrator::Solve(\n+    RightHandSideComputation<Position> compute_acceleration,\n+    SystemState<Position, Variation<Position>> const& initial_value,\n+    Time const& t_final,\n+    Time const& first_time_step,\n+    StepSizeController<Position> step_size_controller,\n+    double const safety_factor,\n+    not_null<Solution<Position, Variation<Position>>*> const solution) const {\n+  using Displacement = Difference<Position>;\n+  using Velocity = Variation<Position>;\n+  using Acceleration = Variation<Velocity>;\n+\n+  // Argument checks.\n+  int const dimension = initial_value.positions.size();\n+  CHECK_EQ(initial_value.momenta.size(), dimension);\n+  CHECK_NE(first_time_step, Time());\n+  Sign const integration_direction = Sign(first_time_step);\n+  if (integration_direction.Positive()) {\n+    // Integrating forward.\n+    CHECK_LT(initial_value.time.value, t_final);\n+  } else {\n+    // Integrating backward.\n+    CHECK_GT(initial_value.time.value, t_final);\n+  }\n+  CHECK_GT(safety_factor, 0);\n+  CHECK_LT(safety_factor, 1);\n+\n+  // Time step.\n+  Time h = first_time_step;\n+  // Current time.\n+  DoublePrecision<Time> t = initial_value.time;\n+\n+  // Position increment (high-order).\n+  std::vector<Displacement> ∆q_hat(dimension);\n+  // Velocity increment (high-order).\n+  std::vector<Velocity> ∆v_hat(dimension);\n+  // Current position.\n+  std::vector<DoublePrecision<Position>> q_hat = initial_value.positions;\n+  // Current velocity.\n+  std::vector<DoublePrecision<Velocity>> v_hat = initial_value.momenta;\n+\n+  // Difference between the low- and high-order approximations of the position.\n+  std::vector<Displacement> q_error_estimate(dimension);\n+  // Difference between the low- and high-order approximations of the velocity.\n+  std::vector<Velocity> v_error_estimate(dimension);\n+\n+  // Current Runge-Kutta-Nyström stage.\n+  std::vector<Position> q_stage(dimension);\n+  // Accelerations at each stage.\n+  // TODO(egg): this is a rectangular container, use something more appropriate.\n+  std::vector<std::vector<Acceleration>> g(stages_);\n+  for (auto& g_stage : g) {\n+    g_stage.resize(dimension);\n+  }\n+\n+  bool at_end = false;\n+  double tolerance_to_error_ratio;\n+\n+  // No step size control on the first step.\n+  goto runge_kutta_nyström_step;\n+\n+  while (!at_end) {\n+    // Compute the next step with decreasing step sizes until the error is\n+    // tolerable.\n+    do {\n+      // Adapt step size.\n+      // TODO(egg): find out whether there's a smarter way to compute that root,\n+      // especially if we make the order compile-time.\n+      h *= safety_factor * std::pow(tolerance_to_error_ratio,\n+                                    1.0 / (lower_order_ + 1));\n+\n+    runge_kutta_nyström_step:\n+      // Termination condition.\n+      Time const time_to_end = (t_final - t.value) - t.error;\n+      at_end = integration_direction * h >= integration_direction * time_to_end;\n+      if (at_end) {\n+        // The chosen step size will overshoot.  Clip it to just reach the end,\n+        // and terminate if the step is accepted.\n+        h = time_to_end;\n+      }\n+\n+      // Runge-Kutta-Nyström iteration; fills |g|.\n+      for (int i = 0; i < stages_; ++i) {\n+        Time const t_stage = t.value + c_[i] * h;\n+        for (int k = 0; k < dimension; ++k) {\n+          Acceleration ∑j_a_ij_g_jk = Acceleration();","path":"integrators/explicit_embedded_runge_kutta_nyström_integrator_body.hpp","position":157,"original_position":157,"commit_id":"cd81b9d28e94d0b85f0c6b89f772fb40b08b8672","original_commit_id":"cd81b9d28e94d0b85f0c6b89f772fb40b08b8672","user":{"login":"eggrobin","id":2284290,"avatar_url":"https://avatars.githubusercontent.com/u/2284290?v=3","gravatar_id":"","url":"https://api.github.com/users/eggrobin","html_url":"https://github.com/eggrobin","followers_url":"https://api.github.com/users/eggrobin/followers","following_url":"https://api.github.com/users/eggrobin/following{/other_user}","gists_url":"https://api.github.com/users/eggrobin/gists{/gist_id}","starred_url":"https://api.github.com/users/eggrobin/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/eggrobin/subscriptions","organizations_url":"https://api.github.com/users/eggrobin/orgs","repos_url":"https://api.github.com/users/eggrobin/repos","events_url":"https://api.github.com/users/eggrobin/events{/privacy}","received_events_url":"https://api.github.com/users/eggrobin/received_events","type":"User","site_admin":false},"body":"Done.","created_at":"2015-05-21T19:44:09Z","updated_at":"2015-05-21T19:44:09Z","html_url":"https://github.com/mockingbirdnest/Principia/pull/531#discussion_r30839415","pull_request_url":"https://api.github.com/repos/mockingbirdnest/Principia/pulls/531","_links":{"self":{"href":"https://api.github.com/repos/mockingbirdnest/Principia/pulls/comments/30839415"},"html":{"href":"https://github.com/mockingbirdnest/Principia/pull/531#discussion_r30839415"},"pull_request":{"href":"https://api.github.com/repos/mockingbirdnest/Principia/pulls/531"}}},"pull_request":{"url":"https://api.github.com/repos/mockingbirdnest/Principia/pulls/531","id":35755122,"html_url":"https://github.com/mockingbirdnest/Principia/pull/531","diff_url":"https://github.com/mockingbirdnest/Principia/pull/531.diff","patch_url":"https://github.com/mockingbirdnest/Principia/pull/531.patch","issue_url":"https://api.github.com/repos/mockingbirdnest/Principia/issues/531","number":531,"state":"open","locked":false,"title":"Explicit ERKN","user":{"login":"eggrobin","id":2284290,"avatar_url":"https://avatars.githubusercontent.com/u/2284290?v=3","gravatar_id":"","url":"https://api.github.com/users/eggrobin","html_url":"https://github.com/eggrobin","followers_url":"https://api.github.com/users/eggrobin/followers","following_url":"https://api.github.com/users/eggrobin/following{/other_user}","gists_url":"https://api.github.com/users/eggrobin/gists{/gist_id}","starred_url":"https://api.github.com/users/eggrobin/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/eggrobin/subscriptions","organizations_url":"https://api.github.com/users/eggrobin/orgs","repos_url":"https://api.github.com/users/eggrobin/repos","events_url":"https://api.github.com/users/eggrobin/events{/privacy}","received_events_url":"https://api.github.com/users/eggrobin/received_events","type":"User","site_admin":false},"body":"Hopefully this integrator is cleaner than the other ones.\r\nI'll reify the ODE / IVP / step size controller in appropriate structs in another CL, for now let's keep the parameters separate.\r\nNote that the integrator appends to the solution without clearing it. This allows for prolongation of a solution vector (by passing the `back` as the initial value).","created_at":"2015-05-19T17:35:01Z","updated_at":"2015-05-21T19:44:09Z","closed_at":null,"merged_at":null,"merge_commit_sha":"e20a0faf670d409b870bfe5979a9f37001a537b2","assignee":null,"milestone":null,"commits_url":"https://api.github.com/repos/mockingbirdnest/Principia/pulls/531/commits","review_comments_url":"https://api.github.com/repos/mockingbirdnest/Principia/pulls/531/comments","review_comment_url":"https://api.github.com/repos/mockingbirdnest/Principia/pulls/comments{/number}","comments_url":"https://api.github.com/repos/mockingbirdnest/Principia/issues/531/comments","statuses_url":"https://api.github.com/repos/mockingbirdnest/Principia/statuses/cd81b9d28e94d0b85f0c6b89f772fb40b08b8672","head":{"label":"eggrobin:explicit-runge-kutta","ref":"explicit-runge-kutta","sha":"cd81b9d28e94d0b85f0c6b89f772fb40b08b8672","user":{"login":"eggrobin","id":2284290,"avatar_url":"https://avatars.githubusercontent.com/u/2284290?v=3","gravatar_id":"","url":"https://api.github.com/users/eggrobin","html_url":"https://github.com/eggrobin","followers_url":"https://api.github.com/users/eggrobin/followers","following_url":"https://api.github.com/users/eggrobin/following{/other_user}","gists_url":"https://api.github.com/users/eggrobin/gists{/gist_id}","starred_url":"https://api.github.com/users/eggrobin/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/eggrobin/subscriptions","organizations_url":"https://api.github.com/users/eggrobin/orgs","repos_url":"https://api.github.com/users/eggrobin/repos","events_url":"https://api.github.com/users/eggrobin/events{/privacy}","received_events_url":"https://api.github.com/users/eggrobin/received_events","type":"User","site_admin":false},"repo":{"id":19360849,"name":"Principia","full_name":"eggrobin/Principia","owner":{"login":"eggrobin","id":2284290,"avatar_url":"https://avatars.githubusercontent.com/u/2284290?v=3","gravatar_id":"","url":"https://api.github.com/users/eggrobin","html_url":"https://github.com/eggrobin","followers_url":"https://api.github.com/users/eggrobin/followers","following_url":"https://api.github.com/users/eggrobin/following{/other_user}","gists_url":"https://api.github.com/users/eggrobin/gists{/gist_id}","starred_url":"https://api.github.com/users/eggrobin/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/eggrobin/subscriptions","organizations_url":"https://api.github.com/users/eggrobin/orgs","repos_url":"https://api.github.com/users/eggrobin/repos","events_url":"https://api.github.com/users/eggrobin/events{/privacy}","received_events_url":"https://api.github.com/users/eggrobin/received_events","type":"User","site_admin":false},"private":false,"html_url":"https://github.com/eggrobin/Principia","description":"N-Body Gravitation and Better Integrators for Kerbal Space Program","fork":true,"url":"https://api.github.com/repos/eggrobin/Principia","forks_url":"https://api.github.com/repos/eggrobin/Principia/forks","keys_url":"https://api.github.com/repos/eggrobin/Principia/keys{/key_id}","collaborators_url":"https://api.github.com/repos/eggrobin/Principia/collaborators{/collaborator}","teams_url":"https://api.github.com/repos/eggrobin/Principia/teams","hooks_url":"https://api.github.com/repos/eggrobin/Principia/hooks","issue_events_url":"https://api.github.com/repos/eggrobin/Principia/issues/events{/number}","events_url":"https://api.github.com/repos/eggrobin/Principia/events","assignees_url":"https://api.github.com/repos/eggrobin/Principia/assignees{/user}","branches_url":"https://api.github.com/repos/eggrobin/Principia/branches{/branch}","tags_url":"https://api.github.com/repos/eggrobin/Principia/tags","blobs_url":"https://api.github.com/repos/eggrobin/Principia/git/blobs{/sha}","git_tags_url":"https://api.github.com/repos/eggrobin/Principia/git/tags{/sha}","git_refs_url":"https://api.github.com/repos/eggrobin/Principia/git/refs{/sha}","trees_url":"https://api.github.com/repos/eggrobin/Principia/git/trees{/sha}","statuses_url":"https://api.github.com/repos/eggrobin/Principia/statuses/{sha}","languages_url":"https://api.github.com/repos/eggrobin/Principia/languages","stargazers_url":"https://api.github.com/repos/eggrobin/Principia/stargazers","contributors_url":"https://api.github.com/repos/eggrobin/Principia/contributors","subscribers_url":"https://api.github.com/repos/eggrobin/Principia/subscribers","subscription_url":"https://api.github.com/repos/eggrobin/Principia/subscription","commits_url":"https://api.github.com/repos/eggrobin/Principia/commits{/sha}","git_commits_url":"https://api.github.com/repos/eggrobin/Principia/git/commits{/sha}","comments_url":"https://api.github.com/repos/eggrobin/Principia/comments{/number}","issue_comment_url":"https://api.github.com/repos/eggrobin/Principia/issues/comments{/number}","contents_url":"https://api.github.com/repos/eggrobin/Principia/contents/{+path}","compare_url":"https://api.github.com/repos/eggrobin/Principia/compare/{base}...{head}","merges_url":"https://api.github.com/repos/eggrobin/Principia/merges","archive_url":"https://api.github.com/repos/eggrobin/Principia/{archive_format}{/ref}","downloads_url":"https://api.github.com/repos/eggrobin/Principia/downloads","issues_url":"https://api.github.com/repos/eggrobin/Principia/issues{/number}","pulls_url":"https://api.github.com/repos/eggrobin/Principia/pulls{/number}","milestones_url":"https://api.github.com/repos/eggrobin/Principia/milestones{/number}","notifications_url":"https://api.github.com/repos/eggrobin/Principia/notifications{?since,all,participating}","labels_url":"https://api.github.com/repos/eggrobin/Principia/labels{/name}","releases_url":"https://api.github.com/repos/eggrobin/Principia/releases{/id}","created_at":"2014-05-01T23:29:34Z","updated_at":"2015-04-07T03:56:10Z","pushed_at":"2015-05-20T19:18:24Z","git_url":"git://github.com/eggrobin/Principia.git","ssh_url":"git@github.com:eggrobin/Principia.git","clone_url":"https://github.com/eggrobin/Principia.git","svn_url":"https://github.com/eggrobin/Principia","homepage":"http://forum.kerbalspaceprogram.com/threads/68502-WIP-Principia-N-Body-Gravitation-and-Better-Integrators-for-Kerbal-Space-Program","size":20260,"stargazers_count":2,"watchers_count":2,"language":"C++","has_issues":false,"has_downloads":true,"has_wiki":true,"has_pages":false,"forks_count":0,"mirror_url":null,"open_issues_count":0,"forks":0,"open_issues":0,"watchers":2,"default_branch":"master"}},"base":{"label":"mockingbirdnest:master","ref":"master","sha":"3d19c2adfe053df58c7df18a4b774600c05dd429","user":{"login":"mockingbirdnest","id":7462691,"avatar_url":"https://avatars.githubusercontent.com/u/7462691?v=3","gravatar_id":"","url":"https://api.github.com/users/mockingbirdnest","html_url":"https://github.com/mockingbirdnest","followers_url":"https://api.github.com/users/mockingbirdnest/followers","following_url":"https://api.github.com/users/mockingbirdnest/following{/other_user}","gists_url":"https://api.github.com/users/mockingbirdnest/gists{/gist_id}","starred_url":"https://api.github.com/users/mockingbirdnest/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/mockingbirdnest/subscriptions","organizations_url":"https://api.github.com/users/mockingbirdnest/orgs","repos_url":"https://api.github.com/users/mockingbirdnest/repos","events_url":"https://api.github.com/users/mockingbirdnest/events{/privacy}","received_events_url":"https://api.github.com/users/mockingbirdnest/received_events","type":"Organization","site_admin":false},"repo":{"id":16634370,"name":"Principia","full_name":"mockingbirdnest/Principia","owner":{"login":"mockingbirdnest","id":7462691,"avatar_url":"https://avatars.githubusercontent.com/u/7462691?v=3","gravatar_id":"","url":"https://api.github.com/users/mockingbirdnest","html_url":"https://github.com/mockingbirdnest","followers_url":"https://api.github.com/users/mockingbirdnest/followers","following_url":"https://api.github.com/users/mockingbirdnest/following{/other_user}","gists_url":"https://api.github.com/users/mockingbirdnest/gists{/gist_id}","starred_url":"https://api.github.com/users/mockingbirdnest/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/mockingbirdnest/subscriptions","organizations_url":"https://api.github.com/users/mockingbirdnest/orgs","repos_url":"https://api.github.com/users/mockingbirdnest/repos","events_url":"https://api.github.com/users/mockingbirdnest/events{/privacy}","received_events_url":"https://api.github.com/users/mockingbirdnest/received_events","type":"Organization","site_admin":false},"private":false,"html_url":"https://github.com/mockingbirdnest/Principia","description":"N-Body Gravitation and Better Integrators for Kerbal Space Program","fork":false,"url":"https://api.github.com/repos/mockingbirdnest/Principia","forks_url":"https://api.github.com/repos/mockingbirdnest/Principia/forks","keys_url":"https://api.github.com/repos/mockingbirdnest/Principia/keys{/key_id}","collaborators_url":"https://api.github.com/repos/mockingbirdnest/Principia/collaborators{/collaborator}","teams_url":"https://api.github.com/repos/mockingbirdnest/Principia/teams","hooks_url":"https://api.github.com/repos/mockingbirdnest/Principia/hooks","issue_events_url":"https://api.github.com/repos/mockingbirdnest/Principia/issues/events{/number}","events_url":"https://api.github.com/repos/mockingbirdnest/Principia/events","assignees_url":"https://api.github.com/repos/mockingbirdnest/Principia/assignees{/user}","branches_url":"https://api.github.com/repos/mockingbirdnest/Principia/branches{/branch}","tags_url":"https://api.github.com/repos/mockingbirdnest/Principia/tags","blobs_url":"https://api.github.com/repos/mockingbirdnest/Principia/git/blobs{/sha}","git_tags_url":"https://api.github.com/repos/mockingbirdnest/Principia/git/tags{/sha}","git_refs_url":"https://api.github.com/repos/mockingbirdnest/Principia/git/refs{/sha}","trees_url":"https://api.github.com/repos/mockingbirdnest/Principia/git/trees{/sha}","statuses_url":"https://api.github.com/repos/mockingbirdnest/Principia/statuses/{sha}","languages_url":"https://api.github.com/repos/mockingbirdnest/Principia/languages","stargazers_url":"https://api.github.com/repos/mockingbirdnest/Principia/stargazers","contributors_url":"https://api.github.com/repos/mockingbirdnest/Principia/contributors","subscribers_url":"https://api.github.com/repos/mockingbirdnest/Principia/subscribers","subscription_url":"https://api.github.com/repos/mockingbirdnest/Principia/subscription","commits_url":"https://api.github.com/repos/mockingbirdnest/Principia/commits{/sha}","git_commits_url":"https://api.github.com/repos/mockingbirdnest/Principia/git/commits{/sha}","comments_url":"https://api.github.com/repos/mockingbirdnest/Principia/comments{/number}","issue_comment_url":"https://api.github.com/repos/mockingbirdnest/Principia/issues/comments{/number}","contents_url":"https://api.github.com/repos/mockingbirdnest/Principia/contents/{+path}","compare_url":"https://api.github.com/repos/mockingbirdnest/Principia/compare/{base}...{head}","merges_url":"https://api.github.com/repos/mockingbirdnest/Principia/merges","archive_url":"https://api.github.com/repos/mockingbirdnest/Principia/{archive_format}{/ref}","downloads_url":"https://api.github.com/repos/mockingbirdnest/Principia/downloads","issues_url":"https://api.github.com/repos/mockingbirdnest/Principia/issues{/number}","pulls_url":"https://api.github.com/repos/mockingbirdnest/Principia/pulls{/number}","milestones_url":"https://api.github.com/repos/mockingbirdnest/Principia/milestones{/number}","notifications_url":"https://api.github.com/repos/mockingbirdnest/Principia/notifications{?since,all,participating}","labels_url":"https://api.github.com/repos/mockingbirdnest/Principia/labels{/name}","releases_url":"https://api.github.com/repos/mockingbirdnest/Principia/releases{/id}","created_at":"2014-02-08T02:55:32Z","updated_at":"2015-05-21T19:05:36Z","pushed_at":"2015-05-20T19:18:25Z","git_url":"git://github.com/mockingbirdnest/Principia.git","ssh_url":"git@github.com:mockingbirdnest/Principia.git","clone_url":"https://github.com/mockingbirdnest/Principia.git","svn_url":"https://github.com/mockingbirdnest/Principia","homepage":"http://forum.kerbalspaceprogram.com/threads/68502-WIP-Principia-N-Body-Gravitation-and-Better-Integrators-for-Kerbal-Space-Program","size":21829,"stargazers_count":74,"watchers_count":74,"language":"C++","has_issues":true,"has_downloads":true,"has_wiki":true,"has_pages":false,"forks_count":7,"mirror_url":null,"open_issues_count":25,"forks":7,"open_issues":25,"watchers":74,"default_branch":"master"}},"_links":{"self":{"href":"https://api.github.com/repos/mockingbirdnest/Principia/pulls/531"},"html":{"href":"https://github.com/mockingbirdnest/Principia/pull/531"},"issue":{"href":"https://api.github.com/repos/mockingbirdnest/Principia/issues/531"},"comments":{"href":"https://api.github.com/repos/mockingbirdnest/Principia/issues/531/comments"},"review_comments":{"href":"https://api.github.com/repos/mockingbirdnest/Principia/pulls/531/comments"},"review_comment":{"href":"https://api.github.com/repos/mockingbirdnest/Principia/pulls/comments{/number}"},"commits":{"href":"https://api.github.com/repos/mockingbirdnest/Principia/pulls/531/commits"},"statuses":{"href":"https://api.github.com/repos/mockingbirdnest/Principia/statuses/cd81b9d28e94d0b85f0c6b89f772fb40b08b8672"}}}}