29 for (Planets::iterator a=tab.begin(); a != tab.end(); a++) {
30 for (Planets::iterator b=tab.begin(); b != a; b++) {
31 double x = a->x - b->x;
32 double y = a->y - b->y;
33 double r = sqrt(x*x + y*y);
34 double f = a->m * b->m * G / (r*r+1);
37 a->vx -= dt * fx / a->m;
38 a->vy -= dt * fy / a->m;
39 b->vx += dt * fx / b->m;
40 b->vy += dt * fy / b->m;
43 for (Planets::iterator a=tab.begin(); a != tab.end(); a++) {
49 Solver(
int n): tab(n), dt(1.0e-4), G(1.0) {
51 for (Planets::iterator a=tab.begin(); a != tab.end(); a++) {
59 for (
int i=0;i<n;i++) Iteration();
70 Rcpp::NumericVector
x,
y,
m,
vx,
vy;
71 for (Solver::Planets::iterator a=s->
tab.begin(); a != s->
tab.end(); a++) {
78 return Rcpp::DataFrame::create(Rcpp::Named(
"x") = x,
80 Rcpp::Named(
"mass") = m,
81 Rcpp::Named(
"Vx") = vx,
82 Rcpp::Named(
"Vy") = vy);
85 if ((
size_t)tab.nrows() != s->
tab.size()) {
88 Rcpp::NumericVector
x = tab[
"x"];
89 Rcpp::NumericVector
y = tab[
"y"];
90 Rcpp::NumericVector
m = tab[
"mass"];
91 Rcpp::NumericVector
vx = tab[
"Vy"];
92 Rcpp::NumericVector
vy = tab[
"Vy"];
93 for (
int i=0;i<tab.nrows();i++) {
110 SEXP
Dollar(Rcpp::XPtr<Wrapper> obj, std::string name) {
111 if (name ==
"data") {
112 return obj->getData();
113 }
else if (name ==
"G") {
114 return Rcpp::NumericVector(1,obj->G());
115 }
else if (name ==
"dt") {
116 return Rcpp::NumericVector(1,obj->dt());
123 Rcpp::XPtr<Wrapper>
DollarAssign(Rcpp::XPtr<Wrapper> obj, std::string name, SEXP v) {
124 if (name ==
"data") {
126 }
else if (name ==
"G") {
127 obj->G() = Rcpp::NumericVector(v)[0];
128 }
else if (name ==
"dt") {
129 obj->dt() = Rcpp::NumericVector(v)[0];
135 Rcpp::CharacterVector
Names(Rcpp::XPtr<Wrapper> obj) {
136 Rcpp::CharacterVector ret;
137 ret.push_back(
"data");
143 int main(
int argc,
char *argv[]) {
146 RInside R(argc, argv,
false,
false,
true);
147 Rcpp::XPtr<Wrapper> wr(
new Wrapper(&S));
148 wr.attr(
"class") =
"Solver";
150 R[
"$.Solver"] = Rcpp::InternalFunction(&
Dollar);
151 R[
"$<-.Solver"] = Rcpp::InternalFunction(&
DollarAssign);
152 R[
"names.Solver"] = Rcpp::InternalFunction(&
Names);
155 std::cout <<
"Want to go interactive? [y/n]";
157 }
while ( !std::cin.fail() && type!=
'y' && type!=
'n' );
159 std::cout <<
"Running an interactive R session. You can explore (and modify) the data in the 'Solver' object." << std::endl;
160 std::cout <<
"[ You can finish the with Ctrl+D ]" << std::endl;
161 R.
parseEval(
"options(prompt = 'R console > ')");
164 std::cout <<
"Finishing the interactive mode" << std::endl;
168 for (
int i=0;i<2000;i++) {
171 R.
parseEval(
"plot(P$x,P$y,cex=P$mass,asp=1,xlim=c(-5,5),ylim=c(-5,5),xlab='X',ylab='Y')");
SEXP Dollar(Rcpp::XPtr< Wrapper > obj, std::string name)
int parseEval(const std::string &line, SEXP &ans)
std::vector< Planet > Planets
Rcpp::CharacterVector Names(Rcpp::XPtr< Wrapper > obj)
Rcpp::DataFrame getData()
int main(int argc, char *argv[])
void setData(Rcpp::DataFrame tab)
Rcpp::XPtr< Wrapper > DollarAssign(Rcpp::XPtr< Wrapper > obj, std::string name, SEXP v)