00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include "orsa_interaction.h"
00026 #include "orsa_secure_math.h"
00027 #include "orsa_universe.h"
00028 #include "orsa_error.h"
00029
00030 #include <cmath>
00031 #include <iostream>
00032
00033 using namespace std;
00034
00035 namespace orsa {
00036
00037 std::string label(const InteractionType it) {
00038
00039 std::string s("");
00040
00041 switch (it) {
00042 case NEWTON: s = "Newton"; break;
00043 case NEWTON_MPI: s = "Newton (MPI)"; break;
00044 case ARMONICOSCILLATOR: s = "Armonic Oscillator"; break;
00045 case GALACTIC_POTENTIAL_ALLEN: s = "Galactic Potential (Allen)"; break;
00046 case GALACTIC_POTENTIAL_ALLEN_PLUS_NEWTON: s = "Galactic Potential (Allen) + Newton"; break;
00047 case JPL_PLANETS_NEWTON: s = "JPL planets + Newton"; break;
00048 case GRAVITATIONALTREE: s = "Gravitational TreeCode"; break;
00049 case RELATIVISTIC: s = "Newton + Relativistic effects"; break;
00050 }
00051
00052 return s;
00053 }
00054
00055 void make_new_interaction(Interaction ** i, const InteractionType type) {
00056
00057 delete (*i);
00058 (*i) = 0;
00059
00060 switch (type) {
00061 case NEWTON: (*i) = new Newton; break;
00062 case NEWTON_MPI:
00063 #ifdef HAVE_MPI
00064 (*i) = new Newton_MPI;
00065 #else
00066 ORSA_WARNING("read NEWTON_MPI interaction from application without MPI support.");
00067 #endif
00068 break;
00069 case ARMONICOSCILLATOR: (*i) = new ArmonicOscillator(1,1); break;
00070 case GALACTIC_POTENTIAL_ALLEN: (*i) = new GalacticPotentialAllen; break;
00071 case GALACTIC_POTENTIAL_ALLEN_PLUS_NEWTON: (*i) = new GalacticPotentialAllenPlusNewton; break;
00072 case JPL_PLANETS_NEWTON: break;
00073 case GRAVITATIONALTREE: (*i) = new GravitationalTree; break;
00074 case RELATIVISTIC: (*i) = new Relativistic; break;
00075 }
00076 }
00077
00078
00079
00080 void MappedTable::load(const std::vector<Body> & f, const bool skip_JPL_planets) {
00081 const unsigned int f_size = f.size();
00082 N = f_size;
00083 mapping.resize(N);
00084 M = 0;
00085 for (unsigned int k=0; k<f_size; ++k) {
00086 mapping[k] = k;
00087 if (!f[k].has_zero_mass()) {
00088 mapping[k] = mapping[M];
00089 mapping[M] = k;
00090 ++M;
00091 }
00092 }
00093
00094 MN = M*N;
00095
00096 if (MN != distance_vector.size()) {
00097 distance_vector.resize(MN);
00098 d1.resize(MN);
00099 d2.resize(MN);
00100 d3.resize(MN);
00101 d4.resize(MN);
00102 one_over_distance.resize(MN);
00103 one_over_distance_square.resize(MN);
00104 one_over_distance_cube.resize(MN);
00105 distance_vector_over_distance_cube.resize(MN);
00106 }
00107
00108 Vector d;
00109 Vector v;
00110 double l;
00111 double one_over_d;
00112 double one_over_d2;
00113 double one_over_d3;
00114 unsigned int index;
00115
00116 for (unsigned int i=0; i<(N-1); ++i) {
00117 for (unsigned int j=i+1; j<N; ++j) {
00118
00119 if (f[i].has_zero_mass() && f[j].has_zero_mass()) continue;
00120 if (skip_JPL_planets && (f[i].JPLPlanet() != NONE) && (f[i].JPLPlanet() != NONE)) continue;
00121
00122 index = ij_to_index(i,j);
00123
00124
00125
00126 d = f[i].DistanceVector(f[j]);
00127
00128 if (distance_vector[index] == d) continue;
00129
00130 l = d.Length();
00131 one_over_d = 1.0/l;
00132 one_over_d2 = one_over_d*one_over_d;
00133 one_over_d3 = one_over_d*one_over_d2;
00134
00135
00136
00137
00138 distance_vector[index] = d;
00139 d1[index] = l;
00140 d2[index] = l*l;
00141 d3[index] = d2[index]*l;
00142 d4[index] = d3[index]*l;
00143 one_over_distance[index] = one_over_d;
00144 one_over_distance_square[index] = one_over_d2;
00145 one_over_distance_cube[index] = one_over_d3;
00146 distance_vector_over_distance_cube[index] = d*one_over_d3;
00147 }
00148 }
00149 }
00150
00151
00152
00153 Newton::Newton() : Interaction(), include_multipole_moments(false), include_relativistic_effects(false), include_fast_relativistic_effects(false), one_over_c2(1.0/(GetC()*GetC())) {
00154 skip_JPL_planets = false;
00155 }
00156
00157 Newton::Newton(const Newton & n) : Interaction(), include_multipole_moments(n.include_multipole_moments), include_relativistic_effects(n.include_relativistic_effects), include_fast_relativistic_effects(n.include_fast_relativistic_effects), one_over_c2(1.0/(GetC()*GetC())) {
00158 skip_JPL_planets = n.skip_JPL_planets;
00159 }
00160
00161 Interaction * Newton::clone() const {
00162 return new Newton(*this);
00163 }
00164
00165 void Newton::fast_newton_acc(const Frame & f, std::vector<Vector> & a) {
00166
00167 const unsigned int size = f.size();
00168
00169 Vector d;
00170
00171 double l;
00172
00173 for (unsigned int i=0;i<size-1;++i) {
00174 for (unsigned int j=i+1;j<size;++j) {
00175
00176 if (zero_mass[i] && zero_mass[j]) continue;
00177
00178 if (skip[i] && skip[j]) continue;
00179
00180 d = f[i].DistanceVector(f[j]);
00181
00182 l = d.Length();
00183
00184 if (d.IsZero()) {
00185 ORSA_WARNING("two objects in the same position! (%s and %s)",f[i].name().c_str(),f[j].name().c_str());
00186 continue;
00187 }
00188
00189 d /= l*l*l;
00190
00191 a[i] += d * f[j].mu();
00192 a[j] -= d * f[i].mu();
00193
00194
00195 }
00196 }
00197 }
00198
00199 void Newton::Acceleration(const Frame & f, vector<Vector> & a) {
00200
00201 const unsigned int size = f.size();
00202
00203 if (size < 2) return;
00204
00205 a.resize(size);
00206
00207
00208
00209
00210
00211 zero_mass.resize(size);
00212 for (unsigned int i=0;i<size;++i) {
00213 zero_mass[i] = f[i].has_zero_mass();
00214 }
00215
00216 skip.resize(size);
00217 if (skip_JPL_planets) {
00218 for (unsigned int k=0;k<size;++k) {
00219 skip[k] = (f[k].JPLPlanet() != NONE);
00220 }
00221 } else {
00222 for (unsigned int k=0;k<size;++k) {
00223 skip[k] = false;
00224 }
00225 }
00226
00227 const bool only_newton = !(include_multipole_moments || include_relativistic_effects || include_fast_relativistic_effects);
00228
00229 if (only_newton) {
00230 for (unsigned int i=0;i<size;++i) {
00231 a[i].Set(0.0,0.0,0.0);
00232 }
00233
00234 fast_newton_acc(f,a);
00235
00236 return;
00237 }
00238
00239 if (include_relativistic_effects && include_fast_relativistic_effects) {
00240 ORSA_WARNING("Both the accurate and the fast version of the Relativistic corrections are activated!");
00241 }
00242
00243 if (a_newton.size() != size) {
00244 a_newton.resize(size);
00245 a_multipoles.resize(size);
00246 a_relativity.resize(size);
00247 }
00248
00249 for (unsigned int i=0;i<size;++i) {
00250 a_newton[i].Set(0.0,0.0,0.0);
00251 a_multipoles[i].Set(0.0,0.0,0.0);
00252 a_relativity[i].Set(0.0,0.0,0.0);
00253 }
00254
00255
00256 Frame mod_f = f;
00257
00258 if (include_relativistic_effects) {
00259
00260
00261
00262
00263
00264
00265 const Vector barycenter = f.Barycenter();
00266 const Vector barycenter_velocity = f.BarycenterVelocity();
00267
00268
00269
00270
00271 for (unsigned int k=0; k<f.size(); ++k) {
00272 mod_f[k].AddToPosition(-barycenter);
00273 mod_f[k].AddToVelocity(-barycenter_velocity);
00274 }
00275
00276 mapped_table.load(mod_f,skip_JPL_planets);
00277 } else {
00278 mapped_table.load(f,skip_JPL_planets);
00279 }
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311 if (include_multipole_moments) {
00312 axis.resize(size);
00313 x_axis.resize(size);
00314 R1.resize(size);
00315 R2.resize(size);
00316 R3.resize(size);
00317 R4.resize(size);
00318
00319 {
00320 Angle alpha,delta,meridian;
00321 Vector planet_axis;
00322 Vector planet_x_axis;
00323 for (unsigned int i=0;i<size;++i) {
00324 if (f[i].JPLPlanet() == NONE) continue;
00325 alpha_delta_meridian(f[i].JPLPlanet(),f,alpha,delta,meridian);
00326
00327 planet_axis.Set(0,0,1);
00328 planet_x_axis.Set(1,0,0);
00329
00330
00331 planet_axis.rotate(meridian.GetRad(),0,0);
00332 planet_x_axis.rotate(meridian.GetRad(),0,0);
00333
00334 planet_axis.rotate(0,halfpi-delta.GetRad(),0);
00335 planet_x_axis.rotate(0,halfpi-delta.GetRad(),0);
00336
00337 planet_axis.rotate(halfpi+alpha.GetRad(),0,0);
00338 planet_x_axis.rotate(halfpi+alpha.GetRad(),0,0);
00339
00340 if (universe->GetReferenceSystem() == ECLIPTIC) {
00341
00342 planet_axis.rotate(0,-obleq_J2000().GetRad(),0);
00343 planet_x_axis.rotate(0,-obleq_J2000().GetRad(),0);
00344 }
00345
00346 axis[i] = planet_axis;
00347 x_axis[i] = planet_x_axis;
00348 }
00349 }
00350
00351 {
00352 double R;
00353 for (unsigned int i=0;i<size;++i) {
00354 R = f[i].radius();
00355 if (R1[i] != R) {
00356 R1[i] = R;
00357 R2[i] = R*R;
00358 R3[i] = R2[i]*R;
00359 R4[i] = R3[i]*R;
00360 }
00361 }
00362 }
00363 }
00364
00365 for (unsigned int i=0;i<size-1;++i) {
00366 for (unsigned int j=i+1;j<size;++j) {
00367
00368 if (zero_mass[i] && zero_mass[j]) continue;
00369
00370 if (skip[i] && skip[j]) continue;
00371
00372 if (mapped_table.Distance(i,j) < (std::numeric_limits<double>::min() * 1.0e3)) {
00373 ORSA_WARNING("two objects in the same position! (%s and %s)",f[i].name().c_str(),f[j].name().c_str());
00374 continue;
00375 }
00376
00377 a_newton[i] += mapped_table.DistanceVectorOverDistanceCube(j,i) * f[j].mu();
00378 a_newton[j] += mapped_table.DistanceVectorOverDistanceCube(i,j) * f[i].mu();
00379 }
00380 }
00381
00382 if (include_relativistic_effects) {
00383
00384 for (unsigned int i=0;i<size-1;++i) {
00385 for (unsigned int j=i+1;j<size;++j) {
00386
00387
00388
00389
00390
00391
00392
00393
00394 if (zero_mass[i] && zero_mass[j]) continue;
00395
00396 if (skip[i] && skip[j]) continue;
00397
00398
00399
00400
00401
00402
00403
00404
00405 double sum_ik = 0.0;
00406 double sum_jk = 0.0;
00407 for (unsigned int k=0; k<size; ++k) {
00408 if (i != k) sum_ik += f[k].mu()*mapped_table.OneOverDistance(i,k);
00409 if (j != k) sum_jk += f[k].mu()*mapped_table.OneOverDistance(j,k);
00410 }
00411
00412 a_relativity[i] += one_over_c2*f[j].mu()*mapped_table.DistanceVectorOverDistanceCube(j,i)*(
00413 - 4.0*sum_ik
00414 - 1.0*sum_jk
00415 + mod_f[i].velocity().LengthSquared()
00416 + 2.0*mod_f[j].velocity().LengthSquared()
00417 - 4.0*mod_f[i].velocity()*mod_f[j].velocity()
00418 - 1.5*pow(mapped_table.DistanceVector(i,j)*mod_f[j].velocity()*mapped_table.OneOverDistance(i,j),2)
00419 + 0.5*mapped_table.DistanceVector(j,i)*a_newton[j]
00420 )
00421 + one_over_c2*f[j].mu()*mapped_table.OneOverDistanceCube(i,j)*((mapped_table.DistanceVector(i,j)*(4.0*mod_f[i].velocity()-3.0*mod_f[j].velocity()))*(mod_f[i].velocity()-mod_f[j].velocity()))
00422 + 3.5*one_over_c2*f[j].mu()*mapped_table.OneOverDistance(i,j)*a_newton[j];
00423
00424 a_relativity[j] += one_over_c2*f[i].mu()*mapped_table.DistanceVectorOverDistanceCube(i,j)*(
00425 - 4.0*sum_jk
00426 - 1.0*sum_ik
00427 + mod_f[j].velocity().LengthSquared()
00428 + 2.0*mod_f[i].velocity().LengthSquared()
00429 - 4.0*mod_f[j].velocity()*mod_f[i].velocity()
00430 - 1.5*pow(mapped_table.DistanceVector(j,i)*mod_f[i].velocity()*mapped_table.OneOverDistance(j,i),2)
00431 + 0.5*mapped_table.DistanceVector(i,j)*a_newton[i]
00432 )
00433 + one_over_c2*f[i].mu()*mapped_table.OneOverDistanceCube(j,i)*((mapped_table.DistanceVector(j,i)*(4.0*mod_f[j].velocity()-3.0*mod_f[i].velocity()))*(mod_f[j].velocity()-mod_f[i].velocity()))
00434 + 3.5*one_over_c2*f[i].mu()*mapped_table.OneOverDistance(j,i)*a_newton[i];
00435 }
00436 }
00437 }
00438
00439 if (include_fast_relativistic_effects) {
00440 for (unsigned int i=0;i<size-1;++i) {
00441 for (unsigned int j=i+1;j<size;++j) {
00442
00443
00444
00445 if (zero_mass[i] && zero_mass[j]) continue;
00446
00447 if (skip[i] && skip[j]) continue;
00448
00449
00450
00451
00452
00453 const Vector r = mapped_table.DistanceVector(i,j);
00454
00455 const Vector v = mod_f[i].velocity() - mod_f[j].velocity();
00456 const double v2 = v.LengthSquared();
00457
00458
00459
00460
00461
00462
00463 a_relativity[i] += one_over_c2*f[j].mu()*mapped_table.OneOverDistanceCube(i,j)*((4.0*f[j].mu()*mapped_table.OneOverDistance(i,j)-v2)*r +
00464 4.0*(r*v)*v);
00465 a_relativity[j] -= one_over_c2*f[i].mu()*mapped_table.OneOverDistanceCube(i,j)*((4.0*f[i].mu()*mapped_table.OneOverDistance(i,j)-v2)*r +
00466 4.0*(r*v)*v);
00467 }
00468 }
00469 }
00470
00471 if (include_multipole_moments) {
00472
00473
00474
00475
00476
00477
00478 for (unsigned int i=0;i<size;++i) {
00479 for (unsigned int j=0;j<size;++j) {
00480
00481 if (i == j) continue;
00482
00483 if (zero_mass[i] && zero_mass[j]) continue;
00484
00485 if (skip[i] && skip[j]) continue;
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496 if ((R1[i] > 0.0) && (fabs(f[i].J2()) > 0.0) && (f[i].JPLPlanet() != NONE)) {
00497
00498
00499
00500 const Vector unit_x_local = mapped_table.DistanceVector(j,i).Normalized();
00501 const Vector unit_y_local = ExternalProduct(axis[i],unit_x_local);
00502 const Vector unit_z_local = ExternalProduct(unit_x_local,unit_y_local);
00503
00504 const Vector tmp_y_axis = ExternalProduct(axis[i],x_axis[i]);
00505 const double tmp_uy_y = unit_y_local*tmp_y_axis;
00506 const double tmp_uy_mx = unit_y_local*x_axis[i];
00507 const double lambda = atan2(tmp_uy_mx,tmp_uy_y);
00508
00509 double s1l,c1l; orsa::sincos(lambda,s1l,c1l);
00510 double s2l,c2l; orsa::sincos(lambda,s2l,c2l);
00511 double s3l,c3l; orsa::sincos(lambda,s3l,c3l);
00512 double s4l,c4l; orsa::sincos(lambda,s4l,c4l);
00513
00514 const double theta = acos(unit_x_local*axis[i]);
00515
00516 double st,ct; orsa::sincos(theta,st,ct);
00517 const double cosec_theta = 1.0/st;
00518
00519 const double ar = R1[i]*mapped_table.OneOverDistance(i,j);
00520 const double ar2 = ar*ar;
00521 const double ar3 = ar*ar2;
00522 const double ar4 = ar*ar3;
00523
00524 const Legendre l(axis[i]*unit_x_local);
00525
00526
00527
00528
00529 const double base_coefficient = - mapped_table.OneOverDistanceSquare(i,j);
00530
00531 const double acc_x_local =
00532 base_coefficient *
00533 (
00534 f[i].J2()*ar2*3.0*l.P2 +
00535 f[i].J3()*ar3*4.0*l.P3 +
00536 f[i].J4()*ar4*5.0*l.P4 +
00537
00538 ar2*(-3.0*l.P22*(f[i].C22()*c2l) ) +
00539 ar3*(-4.0*l.P31*(f[i].C31()*c1l+f[i].S31()*s1l) +
00540 -4.0*l.P32*(f[i].C32()*c2l+f[i].S32()*s2l) +
00541 -4.0*l.P33*(f[i].C33()*c3l+f[i].S33()*s3l) ) +
00542 ar4*(-5.0*l.P41*(f[i].C41()*c1l+f[i].S41()*s1l) +
00543 -5.0*l.P42*(f[i].C42()*c2l+f[i].S42()*s2l) +
00544 -5.0*l.P43*(f[i].C43()*c3l+f[i].S43()*s3l) +
00545 -5.0*l.P44*(f[i].C44()*c4l+f[i].S44()*s4l) ) );
00546
00547 const double acc_y_local =
00548 base_coefficient *
00549 (
00550
00551
00552 ar2*(2.0*cosec_theta*l.P22*(-f[i].C22()*s2l) ) +
00553 ar3*(1.0*cosec_theta*l.P31*(-f[i].C31()*s1l+f[i].S31()*c1l) +
00554 2.0*cosec_theta*l.P32*(-f[i].C32()*s2l+f[i].S32()*c2l) +
00555 3.0*cosec_theta*l.P33*(-f[i].C33()*s3l+f[i].S33()*c3l) ) +
00556 ar4*(1.0*cosec_theta*l.P41*(-f[i].C41()*s1l+f[i].S41()*c1l) +
00557 2.0*cosec_theta*l.P42*(-f[i].C42()*s2l+f[i].S42()*c2l) +
00558 3.0*cosec_theta*l.P43*(-f[i].C43()*s3l+f[i].S43()*c3l) +
00559 4.0*cosec_theta*l.P44*(-f[i].C44()*s4l+f[i].S44()*c4l) ) );
00560
00561 const double acc_z_local =
00562 base_coefficient *
00563 (
00564 f[i].J2()*ar2*(-st)*l.dP2 +
00565 f[i].J3()*ar3*(-st)*l.dP3 +
00566 f[i].J4()*ar4*(-st)*l.dP4 +
00567
00568 ar2*(st*l.dP22*(f[i].C22()*c2l) ) +
00569 ar3*(st*l.dP31*(f[i].C31()*c1l+f[i].S31()*s1l) +
00570 st*l.dP32*(f[i].C32()*c2l+f[i].S32()*s2l) +
00571 st*l.dP33*(f[i].C33()*c3l+f[i].S33()*s3l) ) +
00572 ar4*(st*l.dP41*(f[i].C41()*c1l+f[i].S41()*s1l) +
00573 st*l.dP42*(f[i].C42()*c2l+f[i].S42()*s2l) +
00574 st*l.dP43*(f[i].C43()*c3l+f[i].S43()*s3l) +
00575 st*l.dP44*(f[i].C44()*c4l+f[i].S44()*s4l) ) );
00576
00577 const Vector local_acc =
00578 unit_x_local * acc_x_local +
00579 unit_y_local * acc_y_local +
00580 unit_z_local * acc_z_local;
00581
00582
00583
00584 a_multipoles[i] += f[j].mu() * local_acc;
00585
00586
00587 a_multipoles[j] -= f[i].mu() * local_acc;
00588
00589 }
00590 }
00591 }
00592 }
00593
00594 for (unsigned int i=0;i<size;++i) {
00595 a[i] = a_newton[i] + (a_multipoles[i] + a_relativity[i]);
00596 }
00597
00598
00599 }
00600
00601 double Newton::PotentialEnergy(const Frame &f) {
00602
00603 if (f.size() < 2) return(0.0);
00604
00605 double energy = 0.0;
00606
00607 unsigned int i,j;
00608
00609 Vector d;
00610
00611 double l;
00612
00613 for (i=1;i<f.size();++i) {
00614
00615 if (f[i].mu()==0) continue;
00616
00617 for (j=0;j<i;++j) {
00618
00619 if (f[j].mu()==0) continue;
00620
00621 d = f[i].DistanceVector(f[j]);
00622
00623 l = d.Length();
00624
00625 if (d.IsZero()) {
00626 ORSA_WARNING("two objects in the same position! (%s and %s)",f[i].name().c_str(),f[j].name().c_str());
00627 continue;
00628 }
00629
00630 energy -= f[i].mu()*f[j].mass()/l;
00631
00632 }
00633 }
00634
00635 return (energy);
00636 }
00637
00638
00639
00640 Relativistic::Relativistic() : Interaction(), g(GetG()), c_2(GetC()*GetC()) {
00641
00642 }
00643
00644 Relativistic::Relativistic(const Relativistic &) : Interaction(), g(GetG()), c_2(GetC()*GetC()) {
00645
00646 }
00647
00648 Interaction * Relativistic::clone() const {
00649 return new Relativistic(*this);
00650 }
00651
00652
00653
00654
00655
00656 void Relativistic::Acceleration(const Frame &f, vector<Vector> &a) {
00657
00658 if (f.size() < 2) return;
00659
00660 a.resize(f.size());
00661
00662 unsigned int i,j;
00663
00664 for (i=0;i<a.size();++i)
00665 a[i].Set(0,0,0);
00666
00667 Vector r;
00668 Vector v;
00669
00670 double r_1, r_3;
00671
00672 for (i=1;i<f.size();++i) {
00673
00674 for (j=0;j<i;++j) {
00675
00676
00677
00678 if ((f[i].mass()==0) && (f[j].mass()==0)) continue;
00679
00680
00681 r = f[i].position() - f[j].position();
00682 v = f[i].velocity() - f[j].velocity();
00683
00684 if (r.IsZero()) {
00685 ORSA_WARNING("two objects in the same position! (%s and %s)",f[i].name().c_str(),f[j].name().c_str());
00686 continue;
00687 }
00688
00689 r_1 = r.Length();
00690 r_3 = r.LengthSquared()*r_1;
00691
00692 a[i] -= f[j].mass()/r_3*r + f[j].mass()/(r_3*c_2)*((4.0*g*f[j].mass()/r_1-v.LengthSquared())*r + 4.0*(r*v)*v);
00693
00694
00695
00696 a[j] += f[i].mass()/r_3*r + f[i].mass()/(r_3*c_2)*((4.0*g*f[i].mass()/r_1-v.LengthSquared())*r + 4.0*(r*v)*v);
00697
00698 }
00699
00700 }
00701
00702 for (i=0;i<a.size();++i) a[i] *= g;
00703 }
00704
00705 double Relativistic::PotentialEnergy(const Frame &f) {
00706 Newton newton;
00707 return (newton.PotentialEnergy(f));
00708 }
00709
00710
00711
00712 ArmonicOscillator::ArmonicOscillator(const double free_length_in, const double k_in) : Interaction(), free_length(free_length_in), k(k_in) {
00713
00714 }
00715
00716 ArmonicOscillator::ArmonicOscillator(const ArmonicOscillator & i) : Interaction(), free_length(i.free_length), k(i.k) {
00717
00718 }
00719
00720 Interaction * ArmonicOscillator::clone() const {
00721 return new ArmonicOscillator(*this);
00722 }
00723
00724 void ArmonicOscillator::Acceleration(const Frame& f, vector<Vector>& a) {
00725
00726 if (f.size() < 2) return;
00727
00728 a.resize(f.size());
00729
00730 unsigned int i,j;
00731
00732 Vector d, da;
00733
00734 double ls;
00735
00736 for (i=0;i<f.size();++i)
00737 a[i].Set(0,0,0);
00738
00739 for (i=1;i<f.size();++i) {
00740
00741 if (f[i].mass()==0) continue;
00742
00743 for (j=0;j<i;++j) {
00744
00745
00746
00747
00748
00749
00750
00751
00752 d = f[i].DistanceVector(f[j]);
00753
00754 ls = d.Length();
00755
00756 if (d.IsZero()) {
00757 ORSA_WARNING("two objects in the same position! (%s and %s)",f[i].name().c_str(),f[j].name().c_str());
00758 continue;
00759 }
00760
00761 da = d*(ls-free_length)/ls;
00762
00763 a[i] += da;
00764 a[j] -= da;
00765
00766 }
00767 }
00768
00769 for (i=0;i<a.size();++i)
00770 if (f[i].mass() != 0) a[i] *= k/f[i].mass();
00771
00772 }
00773
00774 double ArmonicOscillator::PotentialEnergy(const Frame &f) {
00775
00776 double energy=0.0;
00777
00778 if (f.size() < 2) return(0.0);
00779
00780 unsigned int i,j;
00781
00782 Vector d;
00783
00784 double ls;
00785
00786 for (i=1;i<f.size();++i) {
00787
00788 if (f[i].mass()==0) continue;
00789
00790 for (j=0;j<i;++j) {
00791
00792
00793
00794
00795
00796
00797 d = f[i].DistanceVector(f[j]);
00798
00799 if (d.IsZero()) {
00800 ORSA_WARNING("two objects in the same position! (%s and %s)",f[i].name().c_str(),f[j].name().c_str());
00801 continue;
00802 }
00803
00804
00805 ls = pow(d.Length()-free_length,2);
00806
00807 energy += ls/2.0;
00808
00809 }
00810 }
00811
00812 return (energy*k);
00813
00814 }
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831 GalacticPotentialAllen::GalacticPotentialAllen() : Interaction() {
00832
00833 g = GetG();
00834
00835 mb = FromUnits(1.40592e10, MSUN);
00836 bb = FromUnits(0.3873e3, PARSEC);
00837 md = FromUnits(8.5608e10, MSUN);
00838 ad = FromUnits(5.3178e3, PARSEC);
00839 bd = FromUnits(0.25e3, PARSEC);
00840 mh = FromUnits(10.3836745e10,MSUN);
00841 ah = FromUnits(12.0e3, PARSEC);
00842
00843 }
00844
00845
00846 GalacticPotentialAllen::GalacticPotentialAllen(const GalacticPotentialAllen &) : Interaction() {
00847
00848 g = GetG();
00849
00850 mb = FromUnits(1.40592e10, MSUN);
00851 bb = FromUnits(0.3873e3, PARSEC);
00852 md = FromUnits(8.5608e10, MSUN);
00853 ad = FromUnits(5.3178e3, PARSEC);
00854 bd = FromUnits(0.25e3, PARSEC);
00855 mh = FromUnits(10.3836745e10,MSUN);
00856 ah = FromUnits(12.0e3, PARSEC);
00857
00858 }
00859
00860 Interaction * GalacticPotentialAllen::clone() const {
00861 return new GalacticPotentialAllen(*this);
00862 }
00863
00864 void GalacticPotentialAllen::Acceleration(const Frame &f, vector<Vector> &a) {
00865
00866 a.resize(f.size());
00867
00868 unsigned int i;
00869 for (i=0;i<a.size();++i) a[i].Set(0,0,0);
00870
00871 double r2,z2,r,rho;
00872 double fbr,fdr,fhr,fhrc;
00873 double fbz,fdz,fhz,fhzc;
00874 double fr,fz;
00875 Vector x;
00876
00877 for (i=0;i<f.size();++i) {
00878
00879
00880
00881 x = f[i].position();
00882
00883
00884 r2 = x.x*x.x+x.y*x.y;
00885 z2 = x.z*x.z;
00886 r = sqrt(r2);
00887 rho = x.Length();
00888
00889
00890 fbr = -(mb*r/(secure_pow(r2+z2+bb*bb,1.5)));
00891 fdr = -(md*r/(secure_pow(r2+secure_pow(ad+sqrt(z2+bd*bd),2),1.5)));
00892 fhr = 1.02*mh*r/(ah*ah*secure_pow(1+1.0/secure_pow(rho/ah,1.02),2)*secure_pow(rho/ah,2.02)*rho);
00893 fhrc = -(mh/(1.02*ah)*(1.0404*(r*secure_pow(rho/ah,0.02))/(ah*rho*secure_pow(1.0+secure_pow(rho/ah,1.02),2))+1.02*(r*secure_pow(rho/ah,0.02))/(ah*rho*(1.0+secure_pow(rho/ah,1.02)))));
00894
00895
00896 fbz = -(mb*x.z/(secure_pow(r2+z2+bb*bb,1.5)));
00897 fdz = -(md*x.z*(ad+sqrt(z2+bd*bd))/(sqrt(z2+bd*bd)*secure_pow(r2+secure_pow(ad+sqrt(z2+bd*bd),2),1.5)));
00898 fhz = 1.02*mh*x.z/(ah*ah*secure_pow(1.0+1.0/secure_pow(rho/ah,1.02),2)*secure_pow(rho/ah,2.02)*rho);
00899 fhzc = -(mh/(1.02*ah)*(1.0404*(x.z*secure_pow(rho/ah,0.02))/(ah*rho*secure_pow(1.0+secure_pow(rho/ah,1.02),2))+1.02*(x.z*secure_pow(rho/ah,0.02))/(ah*rho*(1.0+secure_pow(rho/ah,1.02)))));
00900
00901 fr=fbr+fdr+fhr+fhrc;
00902 fz=fbz+fdz+fhz+fhzc;
00903
00904 a[i].x = fr*x.x/r;
00905 a[i].y = fr*x.y/r;
00906 a[i].z = fz;
00907 }
00908
00909 for (i=0;i<a.size();++i) a[i] *= g;
00910
00911
00912
00913 }
00914
00915 double GalacticPotentialAllen::PotentialEnergy(const Frame &f) {
00916
00917
00918
00919 double energy = 0.0;
00920
00921 double pb,pd,ph;
00922 double r2,z2,r,rho;
00923 unsigned int i;
00924 Vector x;
00925
00926 for (i=0;i<f.size();++i) {
00927
00928 x = f[i].position();
00929
00930
00931 r2 = x.x*x.x+x.y*x.y;
00932 z2 = x.z*x.z;
00933 r = sqrt(r2);
00934 rho = x.Length();
00935
00936 pb = mb/sqrt(r2+z2+bb*bb);
00937 pd = md/sqrt(r2+secure_pow(ad+sqrt(z2+bd*bd),2));
00938 ph = (mh/(rho/ah))*(secure_pow(rho/ah,2.02))/(1.0+secure_pow(rho/ah,1.02));
00939
00940 energy -= (pb+pd+ph);
00941
00942 }
00943
00944 return (energy*g);
00945 }
00946
00947
00948
00949 JPLPlanetsNewton::JPLPlanetsNewton(list<JPL_planets> & l_in) : Interaction(), l(l_in) {
00950 if (universe->GetUniverseType() != Real) {
00951 cerr << "error: using the JPLPlanetsNewton interaction in a non-Real universe!" << endl;
00952 exit(0);
00953 }
00954 g = GetG();
00955 }
00956
00957 JPLPlanetsNewton::JPLPlanetsNewton(const JPLPlanetsNewton & i) : Interaction(), l(i.l) {
00958 if (universe->GetUniverseType() != Real) {
00959 cerr << "error: using the JPLPlanetsNewton interaction in a non-Real universe!" << endl;
00960 exit(0);
00961 }
00962 g = GetG();
00963 }
00964
00965 Interaction * JPLPlanetsNewton::clone() const {
00966 return new JPLPlanetsNewton(*this);
00967 }
00968
00969 void JPLPlanetsNewton::Acceleration(const Frame & f, vector<Vector> & a) {
00970
00971 a.resize(f.size());
00972
00973 if (planets_frame.GetDate() != f.GetDate()) SetupSolarSystem(planets_frame,l,f);
00974
00975 for (unsigned int i=0;i<a.size();++i)
00976 a[i].Set(0,0,0);
00977
00978 for (unsigned int i=0;i<f.size();++i) {
00979 if (f[i].mass() > 0) {
00980 ORSA_ERROR("using the JPLPlanetsNewton interaction with massive objects!");
00981
00982 return;
00983 }
00984 }
00985
00986 Vector d;
00987
00988 double l;
00989
00990 for (unsigned int i=0;i<f.size();++i) {
00991
00992
00993
00994 for (unsigned int j=0;j<planets_frame.size();++j) {
00995
00996
00997
00998
00999 d = f[i].DistanceVector(planets_frame[j]);
01000
01001
01002 l = d.Length();
01003
01004 if (d.IsZero()) {
01005 ORSA_WARNING("two objects in the same position! (%s and %s)",
01006 f[i].name().c_str(),
01007 planets_frame[j].name().c_str());
01008 continue;
01009 }
01010
01011
01012 d /= l*l*l;
01013
01014 a[i] += d * planets_frame[j].mass();
01015
01016 }
01017 }
01018
01019 for (unsigned int i=0;i<a.size();++i) a[i] *= g;
01020 }
01021
01022 double JPLPlanetsNewton::PotentialEnergy(const Frame & f) {
01023 SetupSolarSystem(planets_frame,l,f);
01024 return newton_itg.PotentialEnergy(planets_frame);
01025 }
01026
01027 }