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 #ifndef _ORSA_INTERACTION_H_
00026 #define _ORSA_INTERACTION_H_
00027
00028 #include <vector>
00029 #include <string>
00030
00031 #ifdef HAVE_CONFIG_H
00032 #include <config.h>
00033 #endif // HAVE_CONFIG_H
00034
00035 #include "orsa_coord.h"
00036 #include "orsa_common.h"
00037 #include "orsa_body.h"
00038 #include "orsa_error.h"
00039 #include "orsa_frame.h"
00040
00041 namespace orsa {
00042
00043
00044
00045
00046 enum InteractionType {
00047 NEWTON=1,
00048 ARMONICOSCILLATOR=2,
00049 GALACTIC_POTENTIAL_ALLEN=3,
00050 GALACTIC_POTENTIAL_ALLEN_PLUS_NEWTON=4,
00051 JPL_PLANETS_NEWTON=5,
00052 GRAVITATIONALTREE=6,
00053 NEWTON_MPI=7,
00054 RELATIVISTIC=8
00055 };
00056
00057 inline void convert(InteractionType &it, const unsigned int i) {
00058 switch(i) {
00059 case 1: it = NEWTON; break;
00060 case 2: it = ARMONICOSCILLATOR; break;
00061 case 3: it = GALACTIC_POTENTIAL_ALLEN; break;
00062 case 4: it = GALACTIC_POTENTIAL_ALLEN_PLUS_NEWTON; break;
00063 case 5: it = JPL_PLANETS_NEWTON; break;
00064 case 6: it = GRAVITATIONALTREE; break;
00065 case 7: it = NEWTON_MPI; break;
00066 case 8: it = RELATIVISTIC; break;
00067
00068 default:
00069 ORSA_ERROR("conversion problem: i = %i",i);
00070 break;
00071 }
00072 }
00073
00074 std::string label(const InteractionType it);
00075
00076
00077
00078
00079
00080 class Interaction {
00081 public:
00082 virtual ~Interaction() { };
00083
00084 public:
00085 virtual void Acceleration(const Frame&, std::vector<Vector>&) = 0;
00086 virtual double PotentialEnergy(const Frame&) = 0;
00087
00088 public:
00089 virtual Interaction * clone() const = 0;
00090
00091 public:
00092 virtual bool depends_on_velocity() const { return false; }
00093
00094 public:
00095 void SkipJPLPlanets(const bool b) {
00096 skip_JPL_planets = b;
00097 }
00098 bool IsSkippingJPLPlanets() const {
00099 return skip_JPL_planets;
00100 }
00101 protected:
00102 bool skip_JPL_planets;
00103
00104 public:
00105 virtual InteractionType GetType() const = 0;
00106 };
00107
00108 void make_new_interaction(Interaction**, const InteractionType);
00109
00110
00111
00112 class MappedTable {
00113 public:
00114 void load(const std::vector<Body> & f, const bool skip_JPL_planets);
00115
00116 public:
00117 inline Vector DistanceVector(const unsigned int i, const unsigned int j) const {
00118 return (sign_ij(i,j)*distance_vector[ij_to_index(i,j)]);
00119 }
00120
00121 inline double Distance(const unsigned int i, const unsigned int j) const {
00122 return d1[ij_to_index(i,j)];
00123 }
00124
00125 inline double Distance2(const unsigned int i, const unsigned int j) const {
00126 return d2[ij_to_index(i,j)];
00127 }
00128
00129 inline double Distance3(const unsigned int i, const unsigned int j) const {
00130 return d3[ij_to_index(i,j)];
00131 }
00132
00133 inline double Distance4(const unsigned int i, const unsigned int j) const {
00134 return d4[ij_to_index(i,j)];
00135 }
00136
00137 inline double OneOverDistance(const unsigned int i, const unsigned int j) const {
00138 return one_over_distance[ij_to_index(i,j)];
00139 }
00140
00141 inline double OneOverDistanceSquare(const unsigned int i, const unsigned int j) const {
00142 return one_over_distance_square[ij_to_index(i,j)];
00143 }
00144
00145 inline double OneOverDistanceCube(const unsigned int i, const unsigned int j) const {
00146 return one_over_distance_cube[ij_to_index(i,j)];
00147 }
00148
00149 inline Vector DistanceVectorOverDistanceCube(const unsigned int i, const unsigned int j) const {
00150 return (sign_ij(i,j)*distance_vector_over_distance_cube[ij_to_index(i,j)]);
00151 }
00152
00153 private:
00154
00155
00156
00157 inline unsigned int ij_to_index(const unsigned int i, const unsigned int j) const {
00158 if (swap_ij(i,j)) {
00159 return (mapping[j]*M+mapping[i]);
00160 } else {
00161 return (mapping[i]*M+mapping[j]);
00162 }
00163 }
00164
00165 inline bool swap_ij(const unsigned int i, const unsigned int j) const {
00166 return (i < j);
00167 }
00168
00169 inline double sign_ij(const unsigned int i, const unsigned int j) const {
00170 if (swap_ij(i,j)) {
00171 return -1.0;
00172 } else {
00173 return +1.0;
00174 }
00175 }
00176
00177 inline void index_to_ij(const unsigned int index, unsigned int & i, unsigned int & j) const {
00178 const unsigned int mi = index / M;
00179 const unsigned int mj = index % M;
00180
00181 bool found_i = false;
00182 bool found_j = false;
00183 for (unsigned int k=0; k<N; ++k) {
00184 if (mapping[k] == mi) {
00185 i = k;
00186 found_i = true;
00187 }
00188 if (mapping[k] == mj) {
00189 j = k;
00190 found_j = true;
00191 }
00192 if (found_i && found_j) break;
00193 }
00194 }
00195
00196 private:
00197
00198 unsigned int M, N, MN;
00199
00200 std::vector<unsigned int> mapping;
00201 std::vector<Vector> distance_vector;
00202 std::vector<Vector> distance_vector_over_distance_cube;
00203 std::vector<double> d1;
00204 std::vector<double> d2;
00205 std::vector<double> d3;
00206 std::vector<double> d4;
00207 std::vector<double> one_over_distance;
00208 std::vector<double> one_over_distance_square;
00209 std::vector<double> one_over_distance_cube;
00210 };
00211
00212 class Legendre {
00213 public:
00214 Legendre(const double arg) :
00215 x(arg), x2(x*x), x3(x2*x), x4(x3*x), tmp1(1.0-x2), tmp2(std::sqrt(tmp1)), tmp3(-x/tmp2),
00216 P2(0.5*(3.0*x2-1.0)),dP2(3.0*x),
00217 P3(0.5*(5.0*x3-3.0*x)),dP3(0.5*(15.0*x2-3.0)),
00218 P4(0.125*(35.0*x4-30.0*x2+3.0)),dP4(0.125*(140.0*x3-60.0*x)),
00219 P22(3.0*tmp1),dP22(-6.0*x),
00220 P31(1.5*(1.0-5.0*x2)*tmp2),dP31(1.5*(-10.0*x*tmp2+(1.0-5.0*x2)*tmp3)),
00221 P32(15.0*x*tmp1),dP32(15.0-45.0*x2),
00222 P33(-15.0*tmp1*tmp2),dP33(-15.0*(-2.0*x*tmp2+tmp1*tmp3)),
00223 P41(2.5*x*(3.0-7.0*x2)*tmp2),dP41(2.5*((3.0-21.0*x2)*tmp2+x*(3.0-7.0*x2)*tmp3)),
00224 P42(7.5*(7.0*x2-1.0)*tmp1),dP42(7.5*(14.0*x*tmp1-14.0*x3+2.0*x)),
00225 P43(-105.0*x*tmp1*tmp2),dP43(-105.0*(tmp1*tmp2-2*x2*tmp2+x*tmp1*tmp3)),
00226 P44(105.0*tmp1*tmp1),dP44(-420.0*x*tmp1)
00227 { }
00228
00229 private:
00230 const double x,x2,x3,x4;
00231 private:
00232 const double tmp1, tmp2, tmp3;
00233 public:
00234 const double P2,dP2,P3,dP3,P4,dP4;
00235 const double P22,dP22,P31,dP31,P32,dP32,P33,dP33,P41,dP41,P42,dP42,P43,dP43,P44,dP44;
00236 };
00237
00238
00239
00240 class Newton : public Interaction {
00241 public:
00242 Newton();
00243 Newton(const Newton &);
00244
00245 public:
00246 void Acceleration(const Frame&, std::vector<Vector>&);
00247 double PotentialEnergy(const Frame&);
00248
00249 public:
00250 inline bool depends_on_velocity() const {
00251 return (include_relativistic_effects || include_fast_relativistic_effects);
00252 }
00253
00254 public:
00255 Interaction * clone() const;
00256
00257 InteractionType GetType() const {
00258 return NEWTON;
00259 }
00260
00261 private:
00262 void fast_newton_acc(const Frame &, std::vector<Vector> &);
00263
00264 private:
00265 std::vector<bool> zero_mass;
00266 std::vector<Vector> a_newton;
00267
00268 MappedTable mapped_table;
00269
00270 public:
00271 void IncludeMultipoleMoments(const bool b) {
00272 include_multipole_moments = b;
00273 }
00274 bool IsIncludingMultipoleMoments() const {
00275 return include_multipole_moments;
00276 }
00277 private:
00278 bool include_multipole_moments;
00279
00280 std::vector<Vector> a_multipoles;
00281
00282
00283 std::vector<Vector> axis, x_axis;
00284
00285 std::vector<double> R1;
00286 std::vector<double> R2;
00287 std::vector<double> R3;
00288 std::vector<double> R4;
00289
00290 public:
00291 void IncludeRelativisticEffects(const bool b) {
00292 include_relativistic_effects = b;
00293 }
00294 bool IsIncludingRelativisticEffects() const {
00295 return include_relativistic_effects;
00296 }
00297 private:
00298 bool include_relativistic_effects;
00299 public:
00300 void IncludeFastRelativisticEffects(const bool b) {
00301 include_fast_relativistic_effects = b;
00302 }
00303 bool IsIncludingFastRelativisticEffects() const {
00304 return include_fast_relativistic_effects;
00305 }
00306 private:
00307 bool include_fast_relativistic_effects;
00308
00309 const double one_over_c2;
00310
00311 std::vector<Vector> a_relativity;
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325 private:
00326 std::vector<bool> skip;
00327 };
00328
00329 #ifdef HAVE_MPI
00330 class Newton_MPI : public Interaction {
00331 public:
00332 Newton_MPI();
00333 Newton_MPI(const Newton_MPI &);
00334 ~Newton_MPI();
00335
00336 public:
00337 void Acceleration(const Frame&, std::vector<Vector>&);
00338 double PotentialEnergy(const Frame&);
00339
00340 public:
00341 Interaction * clone() const;
00342
00343 InteractionType GetType() const {
00344 return NEWTON_MPI;
00345 }
00346
00347 private:
00348 double g;
00349
00350 double *bv;
00351 double *av;
00352 double *av_local;
00353 };
00354 #endif // HAVE_MPI
00355
00356 class GravitationalTree : public Interaction {
00357 public:
00358 GravitationalTree();
00359 GravitationalTree(const GravitationalTree &);
00360
00361 public:
00362 void Acceleration(const Frame&, std::vector<Vector>&);
00363 double PotentialEnergy(const Frame&);
00364
00365 public:
00366 Interaction * clone() const;
00367
00368 InteractionType GetType() const {
00369 return GRAVITATIONALTREE;
00370 }
00371
00372 private:
00373 double g;
00374 double theta;
00375 };
00376
00377 class Relativistic : public Interaction {
00378 public:
00379 Relativistic();
00380 Relativistic(const Relativistic &);
00381
00382 public:
00383 void Acceleration(const Frame&, std::vector<Vector>&);
00384 double PotentialEnergy(const Frame&);
00385
00386 public:
00387 Interaction * clone() const;
00388
00389 InteractionType GetType() const {
00390 return RELATIVISTIC;
00391 }
00392
00393 public:
00394 bool depends_on_velocity() const { return true; }
00395
00396 private:
00397 const double g;
00398 const double c_2;
00399 };
00400
00401 class ArmonicOscillator : public Interaction {
00402 public:
00403 ArmonicOscillator(const double free_length_in, const double k_in);
00404 ArmonicOscillator(const ArmonicOscillator &);
00405
00406 public:
00407 void Acceleration(const Frame&, std::vector<Vector>&);
00408 double PotentialEnergy(const Frame&);
00409
00410 public:
00411 Interaction * clone() const;
00412
00413 InteractionType GetType() const {
00414 return ARMONICOSCILLATOR;
00415 }
00416
00417 private:
00418 const double free_length;
00419 const double k;
00420 };
00421
00422 class GalacticPotentialAllen : public Interaction {
00423 public:
00424 GalacticPotentialAllen();
00425 GalacticPotentialAllen(const GalacticPotentialAllen &);
00426
00427 public:
00428 void Acceleration(const Frame&, std::vector<Vector>&);
00429 double PotentialEnergy(const Frame&);
00430
00431 public:
00432 Interaction * clone() const;
00433
00434 InteractionType GetType() const {
00435 return GALACTIC_POTENTIAL_ALLEN;
00436 }
00437
00438 private:
00439 double g;
00440 double mb;
00441 double bb;
00442 double md;
00443 double ad;
00444 double bd;
00445 double mh;
00446 double ah;
00447 };
00448
00449 class GalacticPotentialAllenPlusNewton : public Interaction {
00450 public:
00451 GalacticPotentialAllenPlusNewton() : Interaction() {
00452
00453 }
00454
00455 GalacticPotentialAllenPlusNewton(const GalacticPotentialAllenPlusNewton &) : Interaction() {
00456
00457 }
00458
00459 InteractionType GetType() const {
00460 return GALACTIC_POTENTIAL_ALLEN_PLUS_NEWTON;
00461 }
00462
00463 public:
00464 Interaction * clone() const {
00465 return new GalacticPotentialAllenPlusNewton(*this);
00466 }
00467
00468 inline void Acceleration(const Frame &f, std::vector<Vector> &a) {
00469
00470 tmp_a.resize(a.size());
00471
00472 unsigned int i;
00473
00474 for (i=0;i<a.size();++i)
00475 a[i].Set(0,0,0);
00476
00477 gpa_itg.Acceleration(f,tmp_a);
00478 for (i=0;i<a.size();++i)
00479 a[i] += tmp_a[i];
00480
00481 newton_itg.Acceleration(f,tmp_a);
00482 for (i=0;i<a.size();++i)
00483 a[i] += tmp_a[i];
00484
00485 }
00486
00487 inline double PotentialEnergy(const Frame &f) {
00488 return (gpa_itg.PotentialEnergy(f)+newton_itg.PotentialEnergy(f));
00489 }
00490
00491 private:
00492 GalacticPotentialAllen gpa_itg;
00493 Newton newton_itg;
00494
00495 std::vector<Vector> tmp_a;
00496 };
00497
00498 class JPLPlanetsNewton : public Interaction {
00499 public:
00500 JPLPlanetsNewton(std::list<JPL_planets> &);
00501 JPLPlanetsNewton(const JPLPlanetsNewton &);
00502
00503 public:
00504 void Acceleration(const orsa::Frame &, std::vector<orsa::Vector>&);
00505 double PotentialEnergy(const orsa::Frame &);
00506
00507 public:
00508 Interaction * clone() const;
00509
00510 InteractionType GetType() const {
00511 return JPL_PLANETS_NEWTON;
00512 }
00513
00514 private:
00515 Newton newton_itg;
00516 std::list<JPL_planets> l;
00517 orsa::Frame planets_frame;
00518 double g;
00519 };
00520
00521 }
00522
00523 #endif // _ORSA_INTERACTION_H_