root/sandbox/trunk/pdp/src/N-body.hpp @ 936

/*--------------------------------------------------------------------------
 *  Copyright 2006 Taro L. Saito
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *--------------------------------------------------------------------------*/
//--------------------------------------------------
// N-body.hpp
// Since  2006/10/16 16:06:36
// 
// $Date$ 
// $Author$ 
//--------------------------------------------------
#pragma once
#ifndef __NBODY_HPP20061016161026
#define __NBODY_HPP20061016161026

#include <cmath>
#include <vector>
#include "Triplet.hpp"

/*
 * The class Universe defines global constant values;
 */ 
class Universe
{
public:
        static double gravity() { return double(6.67e-11); }
};


struct Force : public Triplet<double>
{
        Force(double fx, double fy, double fz)  : Triplet<double>(fx, fy, fz) {}

        template <class ForceIter>
        static Force superposition(ForceIter begin, ForceIter end)
        {
                Force f;
                for(ForceIter i = begin; i != end; ++i)
                        f += *i;
                return f;
        }
};

typedef Triplet<double> Velocity;
typedef Triplet<double> Acceleration;


/**
* Particle class represents a two-dimensional point in the universe.
*/
class Particle : public Triplet<double>
{
public:
        Particle(double x, double y, double z, double mass) 
                : Triplet<double>(x, y, z), _mass(mass)
        {}
        virtual ~Particle() {}

        Force gravitationalForce(const Particle& other, double G) const
        {
                Triplet<double> d = *this - other;
                double distSquare = d.x * d.x + d.y * d.y + d.z * d.z;
                double c  = G * _mass * other._mass / (distSquare * std::sqrt(distSquare));
                return Force(c * d.x, c * d.y, c * d.z); // return a projection of F into each of x, y, z-plane.

                // In precise, the above calculation process is as follows:
                //
                // double r = std::sqrt(distSquare);
                // double F = (Universe::gravity() * _mass * other._mass) / (r * r);
                // return Force(F * (dx / r), F * (dy / r), F * (dz / r));      
        }

        template <class ParticleIter>
        void updateAcceleration(ParticleIter begin, ParticleIter end, double G)
        {
                std::vector<Force> force;
                for(ParticleIter i = begin; i != end; ++i)
                {
                        if(this != i)
                                force.push_back(gravitationalForce(*i, G));
                }
                Force currentForce = Force::superposition(force.begin(), force.end());
                _acceleration = currentForce / _mass;   // F = ma
        }

        void move(double dt)
        {
                _velocity += _acceleration * dt;        // V' = V + A * dt
                *this += _velocity * dt;                        // P' = P + V * dt
        }

        double getMass() const { return _mass; } 

        const Velocity& getVelocity() const { return _velocity; }
        void setVelocity(const Velocity& v) { _velocity = v; }
        
        const Acceleration& getAcceleration() const { return _acceleration; }
        void setAcceleration(const Acceleration& a) { _acceleration = a; }
private:
        double _mass;
        Velocity _velocity;
        Acceleration _acceleration;
};


#endif //__NBODY_HPP20061016161026