orsa_interaction.h

Go to the documentation of this file.
00001 /* 
00002    ORSA - Orbit Reconstruction, Simulation and Analysis
00003    Copyright (C) 2002-2004 Pasquale Tricarico
00004    
00005    This program is free software; you can redistribute it and/or
00006    modify it under the terms of the GNU General Public License
00007    as published by the Free Software Foundation; either version 2
00008    of the License, or (at your option) any later version.
00009    
00010    As a special exception, Pasquale Tricarico gives permission to
00011    link this program with Qt commercial edition, and distribute the
00012    resulting executable, without including the source code for the Qt
00013    commercial edition in the source distribution.
00014    
00015    This program is distributed in the hope that it will be useful,
00016    but WITHOUT ANY WARRANTY; without even the implied warranty of
00017    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00018    GNU General Public License for more details.
00019    
00020    You should have received a copy of the GNU General Public License
00021    along with this program; if not, write to the Free Software
00022    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
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   // class Frame;
00044   
00045   // == Unique id == if you need to modify this list, don't change the numbers, but add new numbers
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   // bool depends_on_velocities(const InteractionType it);
00077   
00078   // base class Interaction
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     // i,j are in the frame domain,
00155     // while index is in the mapping domain
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     // M < N, usually is the number of massive particles
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   // derived classes
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     // tmp stuff
00283     std::vector<Vector> axis, x_axis; // axis is the z axis
00284     // std::vector<double> phi; // meridian angle
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        public:
00315        void SkipJPLPlanets(const bool b) {
00316        skip_JPL_planets = b;
00317        }
00318        bool IsSkippingJPLPlanets() const {
00319        return skip_JPL_planets;
00320        }
00321        private:
00322        bool skip_JPL_planets;
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; // 'spring' constant
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; // bulge mass
00441     double bb; // bulge length scale
00442     double md; // disk mass
00443     double ad; // disk length scale
00444     double bd; // disk height scale
00445     double mh; // halo mass (different from the reference article, rescaled to obtain a circular velocity at the sun radius of about 220 km/s)
00446     double ah; // halo length scale
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 } // namespace orsa
00522 
00523 #endif // _ORSA_INTERACTION_H_

Generated on Fri Nov 3 20:37:42 2006 for liborsa by  doxygen 1.4.7